AD-A227  053 


DTIC  FILE  COPY 


(I) 


i  ^ 

^ .  DTIC 

kVv  #%ELECTE| 


SELECTED 
OCT  0  Z 1980 |j 


FREQUENCY-DEPENDENT  MATERIAL  DAMPING 

USING  AUGMENTING  THERMODYNAMIC  FIELDS 

(ATF)  WITH  FRACTIONAL  TIME  DERIVATIVES 

David  S.  Hansen 
First  Lieutenant,  USAF 

AFIT/GAE/ENY/90S-1  ^ 


DEPARTMENT  OF  THE  AIR  FORCE 

AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


pteTHlBUnON  STATEMtJlfrA 

Appr'5''t'.‘l  for  public  ratooMt' 
.  -Mrion  UrUmltud 


Wright-Patlerson  Air  Force  Base,  Ohi( 


90 


099 


y 


(P 


FREQUENCY-DEPENDENT  MATERIAL  DAMPING 
USING  AUGMENTING  THERMODYNAMIC  FIELDS 


(ATF)  WITH  FRACTIONAL  TIME  DERIVATIVES 

David  S.  Hansen 
First  Lieutenant,  USAF 


AFIT/GAE/ENY/90S-1 


DT!C 


ELECTE 
OCT0  21990 


Approved  for  public  release;  distribution  unlimited 


AFIT/GAE/ENY/90S-1 


FREQUENCY -DEPENDENT  MATERIAL  DAMPING  USING 
AUGMENTING  THERMODYNAMIC  FIELDS  (ATF) 
WITH  FRACTIONAL  TIME  DERIVATIVES 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Aeronautical  Engineering 


David  S.  Hansen,  B.S. 
First  Lieutenant,  USAF 


September  1990 


Approved  for  public  release; 


distribution  unlimited 


Preface 


The  purpose  of  this  study  was  to  gain  a  better 
understanding  of  the  modeling  of  material  damping.  With 
increased  operations  in  a  non-atmosphere  environment  such  as 
space,  the  need  exists  for  better  and  more  accurate  models 
of  material  damping,  models  that  include  the  fundamental 
dependence  of  material  damping  on  frequency,  so  that  systems 
designed  for  these  environments  will  be  the  best  possible. 

In  performing  the  research  behind  this  thesis  I 
constantly  relied  on  a  man  who,  in  my  mind,  easily  fits  in 
the  genius  category.  Dr.  Torvik,  I  will  always  be  indebted 


for  your  help,  your  gentle  prodding,  your  knowledge  that 
kept  me  from  floundering  in  deep  waters,  and  your 
friendship.  I  also  wish  to  thank  Lt.  Col.  Bagley,  a  man  I 
greatly  admire,  for  his  enthusiasm  and  overwhelming  desire 
for  my  success,  and  for  his  work  that  helped  this  effort. 

And  to  Lt.  Col.  King  for  his  sage  words  of  advice,  and  help 
as  a  committee  member. 

But  perhaps  I  owe  the  most  to  my  helpmeet  and  eternal 
companion.  Thanks,  Lisa,  for  loving  me  and  pushing  me 
through  the  trying  times  we  had  in  school.  This  degree  is 

truly  yours,  for  you  inspire  me  to  be  better.  Finally,  to  _ 

Fc 

the  rug  rats  who  bring  much  joy  in  my  life.  Nathan,  Ryan, 
and  Megan. . .warm  up  the  Nintendo,  get  the  matchbox  cars  out,  ^ 
practice  up  on  patty-cake  -  Daddy's  coming  home! 
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Abstract 

The  purpose  of  this  study  is  the  development  of  a  model 
using  the  concept  of  augmenting  thermodynamic  fields  (ATF), 
wherein  the  equations  of  thermodynamic  equilibirum  are 
allowed  to  have  derivatives  of  fractional  order.  This 
effort  seeks  to  expand  the  applicabilty  of  the  augmenting 
thermodynamic  field  model,  in  which  augmenting  thermodynamic 
fields  interact  with  the  mechanical  displacement  field.  The 
ATF  model  was  introduced  by  Lesieutre  in  a  dissertation 
submitted  to  the  University  of  California,  Los  Angeles. 

Current  methods  of  analyzing  material  damping  analysis, 
not  being  physically  motivated,  cannot  predict  well  the 
dependency  of  damping  on  frequency.  Two  newer  methods, 
capable  of  predicting  this  frequency  dependence,  are 
discussed.  They  include  the  ATF  model  of  Lesieutre,  and  the 
4-parameter  model,  which  allows  fractional  derivatives  in 
the  description  of  viscoelastic  materials,  introduced  by 
Bagley  and  Torvik  of  the  Air  Force  Institute  of  Technology. 

This  research  effort  applies  the  fractional  derivative 
concept  of  Bagley  and  Torvik' s  model  to  the  model  of 
Lesieutre.  Capitalizing  on  the  previously  made  observation 
that  the  use  of  fractional  derivatives  in  material  models 
enables  the  description  of  material  damping  over  a  much 
broader  frequency  range,  coupled  material  constitutive 


relations  are  developed  using  the  concept  of  augmenting 
thermodynamic  fields,  with  non-integer  differentials  allowed 
in  the  resulting  partial  differential  equations. 

The  complex  modulus  that  results  from  solution  of  these 
partial  differential  equations  is  compared  to  the  complex 
moduli  of  thermoelasticity,  integer-order  viscoelasticity, 
and  viscoelasticity  with  fractional  derivatives  (the 
4-parameter  model)  for  the  case  of  a  uniaxial  rod.  In  each 
case,  the  fractional-order  ATF  model  reduced  to  the 
respective  model,  and,  therefore,  accurately  describes  the 
damping  mechanism  resulting  from  each  of  these  models. 

The  fractional-order  ATF  model  can  be  used  to 
accurately  describe  material  damping  caused  both  by 
thermoelasticity  and  by  viscoelasticity  over  a  broad 
frequency  range. 


FREQUENCY-DEPENDENT  MATERIAL  DAMPING  USING 
AUGMENTING  THERMODYNAMIC  FIELDS  (ATF) 
WITH  FRACTIONAL  TIME  DERIVATIVES 


I . 


Introduction 


Damping  in  Structures 

The  free  vibration  of  any  system  will  eventually  die 
out  due  to  some  form  of  energy  dissipation.  This  energy 
dissipation  is  commonly  referred  to  as  damping.  Damping  in 
a  structure  is  a  result  of  energy  dissipation  mechanisms 
which  can  be  divided  into  three  broad  categories  (Lesieutre, 
1989:4):  (1)  internal  or  material  damping  caused  by  the 
deformation  of  the  material;  (2)  joint  impact  or  friction 


damping,  caused  by  the  dissimilar  deformations  at  structural 
component  interfaces;  and  (3)  energy  dissipation  caused  by 
the  structure's  interaction  with  non-strucutral  components 
or  the  surrounding  environment. 

An  illustration  of  each  of  these  types  of  damping  may 
be  demonstrated  through  an  example  of  a  car  driving  down  the 
highway,  with  the  car  being  the  system  under  consideration. 
Energy  is  dissipated  through  the  movement  and  flexure  of  the 
automobile  components;  this  is  material  damping.  Gears,  the 
engine,  and  interaction  of  the  assembled  components  of  the 
car,  all  represent  joint  or  friction  damping.  Energy  loss 
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caused  by  air  drag  and  the  interaction  of  the  tires  with  the 
highway  are  examples  of  the  third  type  of  damping  of 
structures . 

For  most  systems  in  use  in  the  atmosphere,  the  last  two 
types  of  energy  dissipation  mechanisms  will  dominate.  In 
space,  however,  the  interaction  with  the  environment 
practically  disappears,  and  the  requirements  for  highly 
precise  instruments  and  low  vibration  levels  will  greatly 
reduce  the  loss  mechanism  found  in  joint  and  friction 
damping.  Thus,  material  damping  becomes  much  more 
important,  and  an  accurate  model  of  this  material  damping  is 
needed . 

Material  damping  generally  depends  on  many  factors 
including  frequency,  temperature,  geometry  of  the  structure, 
stress  (or  strain)  history,  and  deformations.  The 
interaction  of  these  factors  leads  to  a  complex  mathematical 
problem.  Zener,  in  his  discussion  of  relaxation  and 
internal  friction  (material  damping),  plots  the  internal 
friction  versus  the  logarithm  of  time,  and  calls  it  the 
relaxtion  spectrum  (Zener,  1948 : Frontispiece) .  Such  a 
figure,  as  replicated  in  figure  1,  is  believed  to  be  typical 
of  the  frequency  dependence  of  material  damping.  Zener 
calculated  the  Internal  friction  caused  by  various  material 
interactions  such  as  grain  boundaries  or  transverse  thermal 
currents,  and  then  linearly  superimposed  the  results  to 
produce  the  graph  of  figure  1.  Many  of  the  damping  models 


Figure  1.  Internal  Friction  (damping)  Versus  Frequency 


in  use  today  are  unable  to  reproduce  this  frequency 
dependence  and  are  therefore  inadequate  in  not  being  fully 
descriptive  of  the  actual  material  behavior. 

Damping  Conventions 

The  world  of  damping  is  not  a  new  one,  although  a  new 
emphasis  has  been  placed  on  it  due  to  expanded  operations  in 
environments  (space)  where  damping  becomes  very  important. 

It  is  important,  therefore,  to  review  the  nomenclature, 
notation,  and  models  that  have  become  standard  in  structural 
dynamics  study. 

Damping  Nomenclature.  The  spring-mass-dashpot  system 
of  figure  2  is  the  simplest  and  most  commonly  used  model  of 
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Figure  2 


Spring-mass-dashpot  System 


systems  with  damping.  It  is  also  known  as  the  viscous 
damping  model,  and  results  in  the  equations  of  motion 

Mx  +  cx  +  Kx  =  F(t)  (1) 

for  the  single  degree-of-freedom  system  with  mass  H, 
stiffness  K,  dashpot  coefficient  c,  and  force  F(t).  Several 
system  parameters  come  from  this  equation.  They  are  listed 
and  briefly  described  below: 

Un  -  Undamped  Natural  Frequency.  Defined  as 
G/n  >  (K/H)  1  2  ,  and  also  known  as  the  natural  frequency,  this 

is  the  frequency  at  which  the  system  would  vibrate  if  there 
were  no  damping  and  no  forcing  function. 

Cc  ~  Critical  Damping.  Defined  as  Cc*2(KH)W2  ,  this 
is  the  dashpot  coefficient  which  would  be  required  for  a 
non-oscil latory  response  of  the  free  system  (F(t)  =  0). 


^  -  Damping  Ratio.  This  is  the  ratio  of  the  actual 
dashpot  coefficent  c  to  the  critical  damping  Cc.  It  can  be 
shown  that  the  time  response  of  a  freely  vibrating 
spring-mass-dashpot  system  displays  a  decay  characterized  by 

e"  S  . 

Damping  Measures.  Several  measures  of  damping  are 
currently  in  use.  Some  are  listed  below,  together  with 
their  definitions. 

Q  -  Resonant  /Amplification  Factor.  The  ratio  of  the 
amplitude  of  the  internal  force  to  that  of  the  external 
force.  The  compliance  is  equal  to  the  reciprocal  of  this. 

At  resonance,  oj  =  ,  and  Q  is  equal  to  one  over  twice 

the  damping  ratio. 

<Jd  “  Logarithmic  Decrement.  This  is  equal  to  the 
natural  logarithm  of  the  ratio  of  the  amplitudes  of  two 
consecutive  displacements  Xi  and  X2,  one  cycle  apart,  or  the 
logarithm  of  the  ratio  of  the  amplitudes  of  two 
displacements  xi  and  xn,  N-cycles  apart  divided  by  the 
number  of  cycles'. 

do  =  (1/N)  ln(xi/XH)  (2) 

This  is  particularly  useful  in  experimentally  finding  the 
damping  of  a  system.  For  small  damping,  d  is  a  small 
quantity,  and  ^  =  d/2 

D  -  Dissipation.  The  dissipation  is  equal  to  the 
amount  of  energy  dissipated  in  one  oscillatory  cycle*. 

-2  rr/cj 

D  =  y  <7  f  dt  ( 3 ) 

0 
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If  divided  by  the  peak  value  of  strain  energy  stored  within 
the  cycle,  Ug,  this  results  in  the  Damping  Capacity,  \J/, 
commonly  used  by  materials  scientists.  The  Damping  Capacity 
is  equal  to  47r  times  the  damping  ratio. 

T]  -  Loss  Factor.  The  ratio  of  the  energy  dissipated 
per  radian  of  vibration  to  the  maximum  energy  stored  is 
equal  to  the  loss  factor.  It  is  also  defined  as  the  ratio 
of  a  material's  loss  modulus  to  its  storage  modulus.  With 
respect  to  the  complex  modulus  E* ,  it  is  the  ratio  of  the 
imaginary  and  real  parts: 

T]  =  I/n(E*  )/J?e(E’ )  (4) 

A  comparison  of  these  damping  measures  to  the  damping 
ratio  results  in  the  following,  for  small  damping: 

1/2Q  =  6:>/27t  =  \J//4n  »  r]/2  (5) 

Damping  Mechanism  Models 

The  linear  matrix  differential  equation  used  most 
commonly  in  structural  dynamics  models  is 

[M]{k}  +  [c]ix}  +  [K]{x}  =  {F}.  (6) 

The  mass  matrix  [M]  is  positive,  definite,  and  the  stiffness 
matrix  [K]  is  positive,  semi-definite.  Both  are  symmetric, 
and  are  easily  found  from  the  mechanics  of  the  problem,  as 
is  the  force  vector  {F} .  The  damping  matrix  [c] ,  although 
positive  semi-definite  and  symmetric,  cannot  be  easily 
derived,  except  in  trivial  cases. 

There  are  several  methods  of  incorporating  damping  into 
structural  models.  These  include,  but  are  not  limited  to. 


viscous  damping,  viscoelastic  damping,  complex  modulus,  and 
structural  damping.  Each  will  be  discussed,  in  scalar  form, 
briefly  below. 

Viscous  Damping.  Viscous  damping  assumes  that  damping 
forces  are  proportional  to  the  time  rate  of  change  of 
strain,  and  thus,  displacement.  This  results  in  the 
equation  of  motion  given  by  equation  (1).  Rewriting  this 
equation  using  the  constants  defined  above  and  setting 
F(t)=0  gives 

x  +  X  +  (Jn  2  X  =  0  (7) 

It  will  later  be  shown  that  the  dissipation  of  a  Kelvin 
solid  increases  monotonical ly  with  frequency.  Most 
materials  do  not  behave  in  such  a  manner,  except  over  small 
frequency  ranges,  and  this  model,  therefore,  although 
mathematically  simple,  is  not  very  useful. 

Viscoelastic  Damping.  The  stress-strain  relationship 
of  a  viscoelastic  material  may  be  described  by  the 
generalized  linear  constitutive  equation  as 

M  d"a  N  d^c 

E  b- -  =  E  an -  (8) 

m=l  dt"  n=l  dtn 

In  the  simplest  case,  this  is  equivalent  to  viscous  damping, 
i  .e. , 

o  -  a.o€  +  aid€/dt  (9) 

since  bo  may  be  set  to  unity  without  loss  of  generality.  In 
this  case,  ao  represents  the  elastic  modulus  ,E,  and  ai 
represents  the  Newtonian  viscosity  coefficient,  . 


The  description  of  more  complex  behavior  is 
accomplished  with  the  inclusion  of  more  terms,  in  turn 
allowing  a  more  accurate  model  over  a  broader  frequency 
range.  This,  however,  quickly  becomes  tedious,  if  not 
impossible  mathematically. 

Complex  Modulus.  Taking  the  Fourier  Transform 

CD 

FT[X(t)]  =  y'e-i‘^tx(t)dt  (10) 

-  CD 

of  equation  (9),  where  temporal  derivatives  introduce  iu), 
gives 

=  E6(tj)  +  i^vU€i(J3)  (11) 

Dividing  both  sides  by  €  gives 

a/e  =  E  +  (12) 

Equation  (12)  suggests  the  possibility  of  a  complex 
modulus,  E*  =  El  +  iEz,  which  may  be  used  in  many  damping 
models.  If  €  =  j?e[eoexp( iwt ) ]  =  €oCos(wt)  for  real  €0, 
then 

=  i?e[E*€  oS  j  ] 

=  f  o  [El cos(a»t )  -  E2Sin(wt)]  (13) 


Recall  that  the  energy  dissipated  per  cycle,  for  real  60,  is 
given  by 


/2  <ii/ 

5 

^  2  (J/ 

j  6o[Ei 


€  dt 


cos(<jt)  -  E2sin(wt) ]  (-w)6oSin(t,;t)dt 


»  E2€o^  TT 


(14) 
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Note  that  D  is  now  independent  of  frequency.  The  peak 
stored  energy.  Us,  is  given  by 

Us  =  1/2  Eito^  (15) 

Therefore,  the  damping  capacity  is  equal  to 

D/Us  =  2  7rE2/Ei  (16) 

and  the  damping  ratio  is 

=  D/(47rU.)  =  1/2  E2/E1.  (17) 

Structural  Damping.  Structural  damping  does  not  depend 
on  the  time  rate  of  strain,  as  in  viscous  damping,  but 
depends  on  the  amplitude  of  oscillation  over  a  wide 
frequency  range.  So,  for  systems  under  harmonic  excitation, 
the  "c"  or  damping  of  equation  (1)  may  be  treated  as 
inversely  proportional  to  the  driving  frequency  (Meirovitch, 
1986:71-73; 1967:400-403) : 

Ceq  =  a/(rr<j)  (18) 

Equation  (1)  may  then  be  rewritten  as 

Mx  +  [a/(.ncj)]x  +  Kx  =  F(t)  (19) 

which  is  not  really  an  equation  since  it  mixes  time 
derivatives  and  frequencies.  The  consequence  is  that  this 
approach  is  only  applicable  in  modeling  constant-amplitude 
forced  vibrations  at  single  frequencies.  Taking  the  Fourier 
transform  of  x(t)  gives  iwX,  where  X  is  the  Fourier 
transform  of  x(t),  thus  motivating  the  rewriting  of  equation 
(19)  as  the  pseudo-equation 

Mx  +  K(1  +  iy)x  =  F(t)  (20) 

where 


y  «  a/(7TK) 
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(21) 


which  is  equivalent  to  the  loss  factor.  The  quantity 
K(1  +  iy)  is  sometimes  called  the  complex  stiffness  or 
complex  damping  model. 

Augmenting  Thermodynamic  Field  model  with  Fractional 
Derivatives .  This  research  effort  applies  concepts  of 
recent  research  undertaken  by  Torvik  and  Bag ley  of  the  Air 
Force  Institute  of  Technology,  Wright  Patterson  Air  Force 
Base,  to  the  research  efforts  of  Lesieutre,  formerly  of  the 
University  of  California,  Los  Angeles.  Bagley  and  Torvik 
have  described,  in  a  series  of  papers  (Bagley  and  Torvik, 
1986;  1983;  Torvik  and  Bagley,  1987;  1984),  the  application 
of  material  models  involving  the  fractional  derivative  to 
problems  in  vibration.  Lesieutre  in  his  doctoral 
dissertation  (Lesieutre,  1989)  introduced  the  notion  of 
augmenting  thermodynamic  fields  to  describe  a  new  method  of 
modeling  frequency-dependent  material  damping  in  structural 
damping  analysis.  This  effort  will  apply  fractional 
derivatives  to  the  augmenting  thermodynamic  field  model  in 
an  effort  to  expand  the  applicability  of  that  model. 

The  utility  of  the  ATF  model  with  fractional 
derivatives  will  be  demonstrated  through  comparison  of  the 
model  with  other  material  damping  models  including 
thermoelasticity,  viscoelasticity,  and  viscoelasticity  with 
fractional  derivatives.  It  will  be  seen  that  the  ATF  model 
with  fractional  derivatives  can  be  reduced  to  these  models, 
and  is,  therefore,  a  more  powerful  model  for  use  in  material 
damping  study. 
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II.  Augmenting  Thermodynamic  Fields 


Noting  the  physical  significance  of  the  internal  state 
variables  of  materials  science,  Lesieutre  (Lesieutre,  1989) 
introduced  augmenting  thermodynamic  fields  to  interact  with 
the  mechanical  displacement  field  of  continuum  dynamics.  He 
then  used  nonequilibrium,  irreversible  thermodynamics  to 
develop  coupled  constitutive  and  partial  differential 
equations.  His  approach  will  be  reviewed  here  as  a  basis 
for  the  work  of  this  effort. 

Anelasticity  (Nowick  and  Berry,  1972:1-5).  Nowick  and 
Berry  intoduce  anelastic  behavior  formally  by  first 
considering  an  ideal  elastic  material,  for  which  the 
relation  of  stress  to  strain  is  defined  by  Hooke's  law 
<7  =  Ec  ,  where  E  is  called  the  modulus  of  elasticity.  For 
this  material,  three  conditions  result.  They  are:  (1)  the 
applied  stress  produces  a  unique  strain  response  through  an 
equilibrium  relationship,  i.e.,  there  is  complete 
recoverability,  (2)  the  equilibrium  relationship  produces  an 
instantaneous  response,  and  (3)  the  response  is  linear.  By 
lifting  the  instantaneity  condition,  the  type  of  behavior 
known  as  anelasticity  is  produced.  With  an  anelastic 
material,  a  time  dependent  nonelastic  response  will  occur  in 
addition  to  the  instantaneous  elastic  reaction. 

All  materials  that  qualify  as  thermodynamic  solids 
satisfy  the  first  condition  above,  since  these  solids  can 
continuously  assume  unique  equilibrium  states  in  response  to 
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small  changes  In  some  external  variable  such  as  stress  or 
temperature.  Relaxation  is  the  self-adjustment  with  time  of 
a  thermodynamic  system  toward  a  new  equilibrium  values  as  a 
result  of  a  change  in  some  external  variable  (Nowick  and 
Berry,  1972:4).  Anelastic  relaxation  is  the  tendency  of  a 
thermodynamic  system  to  evolve  towards  an  equilibrium  state 
when  that  state  depends  on  an  external  mechanical  value  such 
as  stress  or  strain. 

Nowick  and  Berry  summarize  anelastic  relaxation  as  the 
thermodynamic  phenomenon  which  results  from  the  changing  of 
internal  variables  which  couple  stress  and  strain  to  new 
equilbrium  values  through  such  processes  as  kinetic 
diffusion  (Nowick  and  Berry,  1972:5).  Material  damping  is  a 
direct  result  of  such  anelastic  relaxation  processes  in  real 
materials . 

When  the  nature  of  the  associated  thermodynamic 
relaxation  processes  is  considered,  two  classes  of  material 
damping  mechanisms  can  be  defined.  They  are  damping  as  a 
result  of  system  non-uniformity  and  damping  as  a  result  of  a 
local  departure  from  thermodynamic  equilibrium  (Lesieutre, 
1989:23-26) . 

System  Non-Uniformity.  Thermoelasticity  is  a 
we 1 1 -understood  example  of  this  type  of  damping  mechanism. 
Torvlk  reviews  this  mechanism,  stating  that  thermoelastic 
damping  is  a  result  of  the  coupling  created  when  a  thermal 
expansion  term  appears  in  the  equation  of  state  of  a 
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material,  as  well  as  in  the  energy  equation  which  describes 
the  process  (Torvik,  1989:2). 

Consider  an  axial  rod  subjected  to  longitudinal  elastic 
waves.  During  a  cycle,  one-half  each  wavelength  of  the  rod 
undergoes  an  expansion  with  an  accompanying  increase  in 
strain,  and  thus,  an  increase  in  volume.  To  satisfy  the 
energy  balance,  temperature  of  this  half  of  the  rod  must 
decrease.  In  an  opposite  sense,  the  temperature  of  the 
other  half  of  the  rod  is  Increasing,  due  to  the  decrease  in 
volume  (compression)  and  strain  of  that  half  of  the  rod. 

The  end  result  is  a  flow  of  heat  by  conduction  along  the 
longitudinal  axi^:,  causing  temperature  changes  with  time 
along  the  rod.  These  changes  couple  into  the  strain  of  the 
constituitive  equation  causing  material  damping,  although, 
as  was  noted  in  Torvik's  paper,  this  damping  is  small. 

Local  Departure  from  Thermodynamic  Equilibrium 
(Lesieutre,  1989:24).  Also  known  as  microstructural 
damping,  this  damping  is  a  result  of  the  coupling  of  stress 
and  strain  with  characteristics  of  real  materials  other  than 
temperature.  Material  damping  is  affected  by  thermodynamic 
relaxations  due  to  point  defects,  dislocations,  or  grain 
boundaries,  which  are  thermally-activated,  and  are  sensitive 
to  the  composition  and  microstructure  of  the  material . 

These  are  all  internal  variables,  evolving  jointly  with 
stress  and  strain,  which  motivated  the  introduction  of 
augmenting  thermodynamic  fields  (ATF)  by  Lesieutre.  These 
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ATFs  are  to  interact  with  the  usual  mechanical  displacement 
field. 

Irreversible  Thermodynamics.  The  point  defects  of 
microstructural  damping  are  thermally  activated,  leading  one 
to  the  conclusion  that  material  damping  is  fundamentally  a 
thermodynamic  phenomenon.  Kovalenko  states  that  the 
irreversible  thermoelastic  deformation  of  a  body  is  a 
non-equilibrium  process,  with  the  irreversibl ity  being  due 
to  the  temperature  gradient  (Kovalenko,  1969:7).  A  process 
associated  with  an  increase  in  the  entropy,  with  all  real 
processes  falling  within  this  category,  is  said  to  be 
irreversible.  This  increase  in  entropy  provides  a 
quantitative  measure  of  the  irreversibl ity  of  the  process 
(Callen,  1960:63-65). 

In  irreversible  thermody*.amics  the  entropy  balance 
equation  plays  a  central  role  (Lesieutre,  1989:25-27).  The 
rate  of  local  entropy  production  equals  the  entropy  leaving 
the  region,  plus  the  rate  of  entropy  produced  within  the 
region  (Callen,  1960:204-288).  In  others  words,  entropy  of 
a  volume  element  changes  with  time  due  to  (1)  entropy  flow 
into  the  volume  element,  and  (2)  entropy  created  due  to  the 
irreversible  processes  inside  the  element.  The  entropy 
created,  or  the  rate  of  production  of  entropy,  is  the  sum  of 
products  of  each  flux  with  its  associated  affinity.  Fluxes 
characterize  some  irreversible  relaxation  process,  such  as 
in  a  heat  flux.  The  affinity  can  be  thought  of  as  the 


generalized  force  which  drives  the  irreversible  process 
toward  the  equilibrium  state. 

Axial  Rod  With  One  Augmenting  Thermodynamic  Field  (ATF) 
(Lesieutre,  1989:29-41).  Lesieutre  begins  by  considering 
the  case  of  one-dimensional  motion,  corresponding  to 
longitudinal  vibration  of  a  slender  rod  .  Displacement 
along  the  rod  is  denoted  by  u  =  u(x,t)  ,  mass  density, 
which  is  uniform,  by  p,  and  the  unrelaxed  modulus  of 
elasticity  by  E.  He  then  introduces  a  single  augmenting 
thermodynamic  field  (ATF),  ^  =  ^(x,t)  .  The  stress  field, 

a,  is  thermodynamically  conjugate  to  the  strain,  e,  as  is 
the  affinity  field.  A,  to  the  ATF,  From  irreversible 

thermodynamics,  the  affinity  can  be  interpreted  as  the 
generalized  force  which  drives  ^  to  equilibrium.  The 
strength  of  the  coupling  of  the  two  dependent  fields,  u  and 
is  described  by  the  material  property  (5.  If  there  were 
no  coupling,  the  stress  would  be  related  to  strain  by 
Hooke's  law.  Analogously,  the  material  property  a  relates 
changes  in  A  to  those  in 

Helmholtz  Free  Energy.  Spanner  states  that,  to  have 
spontaneity  of  a  process  under  isothermal  conditions,  the 
process  must  involve  a  decrease  in  free  energy  (Spanner, 
1964:Ch  7).  The  Helmholtz  free  energy  is  used  when  no 
external  volume  work  occurs  as  the  process  takes  place. 

The  free  energy  of  a  system  is  not  a  part  of  its  total 
energy.  In  truth,  it  is  not  simply  an  energy  at  all,  but 
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represents  the  potentiality  of  the  system  doing  work  under 
the  condition  that  the  system  is  supplied  with  all  the 
heating  or  cooling  it  requires  to  keep  a  constant 
temperature.  When  the  system  actually  does  work,  it  is  at 
the  expense  of  its  free  energy,  with  part  of  that  work  done 
out  of  its  own  resources,  and  part  out  of  the  heat  energy 
supplied  to  it  from  its  surroundings,  resulting  in  an 
increase  in  system  entropy  (Spanner,  1964:91).  Callen 
states  that  conditions  are  such  in  many  processes  that  the 
ambient  atmosphere  acts  as  a  heat  reservoir  to  maintain  the 
temperature  constant  (Callen,  1960:107).  It  is  for  these 
processes  that  the  Helmholtz  potential  representation  is 
well  suited. 

The  state  of  system  equilibrium  is  determined  in  an 
isothermal  system  by  the  condition  of  minimum  free  energy, 
di .  Here,  di  represents  the  available  useful  work,  other 
than  pressure-volume  work,  at  constant  temperature  and 
volume.  If  df  is  negative  the  system  will  do  work 
spontaneously.  Processes  for  which  this  is  the  case  are 
natural  processes  in  the  sense  that  work  is  done  by  the 
system  in  an  effort  to  approach  equilibrium  (Penner, 

1968:71) . 

The  Helmholtz  free  energy  density,  f,  being  a 
thermodynamic  potential  appropriate  for  use  when  strain  is 
an  independent  variable  (Lesieutre,  1989:31),  can  be  used  to 


derive  the  constitutive  equations: 

f  =  l/2Ee2  i/2a^2  (22) 

Material  Constitutive  Equations.  The  constitutive 
equations  are  then  found  from  the  thermodynamic  relations 
(Spanner,  1964:162) 

o  =  dlldi  ,  A  =  -  dl/d^  (23) 

which  result  in 


a  =  E€  -  ^  ^  (24) 

A  =  <56  -  (25) 

In  classical  irreversible  thermodynamics,  it  is  assumed 
that  the  time  rate  of  change  of  4  is  proportional  to  A,  the 
affinity  or  thermodynamic  force.  For  the  quadratic 
dependence  of  equation  (22)  it  is  proportional  to  its 
deviation  from  an  equilibrium  value,  with  the 
proportionality  constant,  B,  being  a  material  property 
different  for  any  distinct  dissipative  mechanism: 

\  =  -B(^  -  |)  (26) 

where  the  dot  (■)  denotes  temporal  derivatives. 

The  value  of  \  is  the  value  of  ^  at  equilibrium,  or 
when  the  affinity,  the  "generalized  force"  which  drives  ^  to 
equilibrium,  is  zero.  Setting  equation  (25)  to  zero  and 
solving  for  ^  gives 

i  =  bi/oL  (27) 

which  in  turn  can  be  substituted  into  equation  (26)  to  give 

\  =  -B(^  -  (5€/a)  (28) 
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Cauchy's  first  law  is  invoked  in  order  to  insure 
conservation  of  linear  momentum.  In  general 

tik.i  +  p(fk  ~  ak)  =  0  (29) 

where  fk  and  ak  are  components  of  the  body  force  and 
acceleration,  repectively.  For  the  case  of  axial  motion  of 
a  one-dimensional  bar,  and  neglecting  body  forces,  this 
reduces  to 

do/d'x.  -  pa  =  da/d'x.  -  pu  =  0  (30) 

where  the  dot  denotes  temporal  derivatives. 

Recall  that  f(x,t)  =c>u(x,t)/0x  =  u'  ,  where  the  prime 
.denotes  spatial  derivatives.  Using  this,  substituting 
equation  (24)  into  equation  (30),  and  rearranging  terms 
gives,  together  with  equation  (28),  a  set  of  two  coupled 
partial  differential  equations  in  u  and 
pii  -  Eu"  =  -6^' 

^  =  (B(5/a)u'  (31) 

Lesieutre  uses  the  principles  of  nonequilibrium 
irreversible  thermodynamics  to  determine  the  conditions 
under  which  a  given  problem  involving  these  equations  is 
well-posed  (Lesieutre,  1989:32-34).  These  conditions  are 
P  >  0,  E  >  0,  and  a  >  ((Jz/E)  >  0  (32) 

conditions  easily  satisfied  in  practice. 

Four i er  Ana lysis.  The  set  of  equations  (31)  is  a  set 
of  linear  partial  differential  equations,  and,  as  such, 
Fourier  analysis  may  be  used  for  their  solution,  especially 
since  the  frequency-dependence  of  the  damping  is  of  primary 


interest.  Consider  a  plane  wave  traveling  in  the  negative 
x-direction  with  temporal  frequency  u;  and  spatial  frequency 
(wave  number)  k  according  to  ei<*'**wt).  The  Fourier 
transform  of  spatial  derivatives  will  then  introduce  ik,  and 
the  transform  of  time  derivatives  will  introduce  io), 
boundary  conditions  not  being  considered.  Equation  (31) 
yields 

(-pu;2  +  k2E)U  +  ik(5Z  =  0  (33) 

(Bik(5/a)U  -  (B  +  iw)Z  =  0  (34) 

where  U  and  Z  are  the  Fourier  transforms  of  u(x)  and  ^(x), 
respectively . 

Damping  Ratio.  Solving  equation  (34)  for  Z,  and 
substituting  into  equation  (33)  gives,  after  algebraic 
manipulation 

r  1  -  d  +  i(w/B)' 

-  =  E  -  (35) 

k2  L  1  +  i(<^/B) 

where 

A  =  dVEa  (36) 

The  term  inside  the  [  ]  brackets  can  be  interpreted  as 
a  complex  modulus,  E’'(u)).  which  after  further  manipulation 
becomes 

1  -  d  +  (w/B)2  +  i(du>/B) 

E*(a;)  «  E  -  »  0  (37) 

1  +  (w/B)? 

and  the  loss  factor,  rj ,  is  found  by 

rj  m  Iin(E*(ci>)]/l?e(E’(u>)]  (38) 
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or,  from  equation  (37) 

ri  «4(u)/B)/[l  -  4  +  (u>/B)2]  (39) 

If  small  (lamping  is  assumed,  the  damping  ratio  equals 
half  of  the  loss  factor  (equation  (5)): 

^  «  d(w/B)/{2tl  -  d  +  (u)/B)2]>  (40) 

For  small  damping  the  damping  ratio  is  small,  and  the 
numerator  of  equation  (40)  must  be  small  in  comparison  to 
its  denominator.  This  being  the  case,  d  can  be  assumed 
small,  and  can  be  cancelled  from  the  denominator,  leaving 
the  following  ecjuation  for  the  damping  ratio: 

^  »  ((5  V4Ea)(2u)/B)/(l  +  (ai/B)2  1  (41) 

Figure  3  is  a  plot  of  the  loss  factor  versus  frequency 
for  material  property  values  that  correspond  roughly  to 
aluminum  (Lesieutre,  1989:4B): 


Figure  3.  Loss  Factor  vs.  Frequency  for  Single  ATF  Model 


E  =  7.13  X  1010  N/m? 

B  =  8000  1/s 

a  =  8000  N/m2 

6  =  4.7766  X  10<>  M/m2 


As  can  be  seen,  damping  peaks  at  a  frequency  of  (j  =  B  and 
tends  toward  zero  at  both  high  and  low  frequencies. 

The  power  of  the  ATF  model  in  the  description  of 
material  damping  can  be  seen  as  it  is  compared  to  other 
material  damping  models.  It  will  be  seen  that  the  ATF  model 
with  no  fractional  derivatives  can  describe  damping  as 
modeled  by  the  thermoelastic  and  classic  viscoelastic 
damping  models. 


Comparison  to  the  Thermoelastic  Model.  The  linearized 
equations  of  dynamic  coupled  thermoelasticity  are  given  in 
three  dimensions  as  (Boley  and  Weiner,  1960:38) 

<7ii  “OfT(3A'^  2  fj)  6  i  i  {T  ~  To)  (42) 

kT,ii  =  CvT  +  (3A  +  2/;)aTTo€kk  (43) 

with  equilibrium  requiring 

<7  i  i  ,  i  =  p  u  i  ( 44 ) 

and  compatibility  conditions  giving 

€ii~l/2(Ui,i+Ui,i)  (45) 

where 


<7i  i  =  stress 
€ i i  =  strain 
<5  i  i  =  Kronecker  delta 
p  =  mass  density 

T  =  change  in  temperature  from  the  initial , 
stress-free  state 

To  =  initial  temperature,  a  constant 
Ui  =  displacement 
f<  =  thermal  coductivity 
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Cv  =  specific  heat  at  constant  volume  (per  unit 
volume) 

X,  U  =  isothermal  Lame's  constants 

otr  =  coefficient  of  thermal  expansion 

Given  a  one-dimensional  displacement  with  no  shear  the 
following  conditions  result: 

Ui  =  Ui (x, t) 
u?=  U3  =  0 

0^2  =<^33=0  (46 ) 

Compatibility,  equation  (45),  results  in 

6ll  =  du.\/dx  €77-^33  -  0  (47) 

and  equilibrium  gives 

pvi  =  =  dOx/d'x.  (48) 

Following  application  of  the  conditions  of  equation 
(46)  to  equation  (42)  we  have 

<7i  1  =  (A  +  2p)€i  1  -  arOA  +  2p)(T  -  To)  (49) 

Equilibrium,  equation  (44),  requires  the  derivative  of 
equation  (49)  with  respect  to  x.  This,  together  with  the 
result  of  applying  equation  (46)  to  equation  (43),  results 
in  the  linearized  equations  of  dynamic  coupled 
thermoelasticy  for  plane  strain,  as  given  by  Boley  and 
Weiner  (Boley  and  Weiner,  1960:45-73): 

kT"  -  CvT  -  (3A  +  2p)aTToU’  =  0 
(A  +  2p)u"  -  pvi  -  (3A  +  2fi)arT'  =  0  (50) 

Note  that  boundary  condition  effects  are  not  included  in 
this  analysis. 

For  the  case  of  the  uniaxial  rod  considered  in  this 
thesis  a  set  of  conditions  different  than  those  of  equations 
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(46)  and  (47)  result: 

(7??  =  =  0 

Ui  /  0,  U?  =  0,  =  0 

€  1  ]  =  dU]  /d'X. 

€??=e3i=“i^€ii  (51) 

where  i/  is  Poisson's  ratio.  Equation  (42)  then  gives 

dll  =  A  ( ^  1  1  +  ^  ?  2  +  ^3  3 )  +  2H(  \  \  -  01 T  ( T-To  )  (  3A  +  2// )  (52) 

Applying  the  last  of  the  conditions  of  equation  (51)  results 
in  the  equation  of  stress  in  the  x-direction: 

dll  =  A€,](1  -  2^/)  +  2//0i,  -0<t(T  -  To)(3A+  2^)  (53) 

Likewise,  equation  (43)  results  in 

KT"  =  CvT  +  (3A  +  2^)a^Tr,(l  -  2u)\l'  (54) 

Again,  equilibrium  is  required,  which  results  in  the 
following  from  equation  (53): 

(k  +  2^  -  2t/)u"  -  aT(3A  +  2//)T’  -  pii  =  0  (55) 

Recalling  the  definition  of  the  isothermal  Lame's 
constants,  which  involve  the  isothermal  modulus,  Et , 

A  =  (i/Et)/[(1  +  i')  (1  -  2t/) 

p  =  Et/[2(1  +  i^)]  (56) 

and  substitution  of  equation  (56)  into  equations  (53)  and 
(55)  results  in  the  equations  of  dynamic  coupled 
thermoelasticity  for  the  plane  stress  case  of  the  uniaxial 
rod : 

Kl"  -  Cvt  -  EraxToU'  =  0 

Eju"  -  pu  -  EjaxT'  =  0  (57) 
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Taking  the  Fourier  transform  of  equation  (57),  with 
spatial  derivatives  introducing  ik  and  temporal  derivatives 


introducing  iu  ,  results  in  the  following  equations: 

A<-(ik)2T  -  Cv(iw)T  -  ETOerToCik)  (ia>)U  =  0  (58) 

ET(ik)2U  -  p(i<j)2U  -  ETOtT(ik)T  =  0  (59) 

where  U  is  the  Fourier  transform  of  u(x,t)  and  f  is  the 
Fourier  transform  of  T(x,t).  Solving  equation  (59)  for  T, 
and  substituting  the  result  into  equation  (58)  gives,  after 
some  algebraic  manipulation. 


pcj'^ 

k2 


=  El 


(iwToEiaT^  )/(/^k2 ) 


1  + 


(60) 


1  +  ( icjc  V ) /(^k2 ) 

This  is  in  the  form  of  the  familiar  complex  modulus,  with 
the  right  hand  side  of  equation  (60)  being  that  modulus, 

1  +  i[(j(cv  +  aT2ToET)/(Kk2)l 


E* (u)  =  Et 


1  +  ici;(cv/(^k2 )  ] 


(61) 


which  in  turn  can  be  compared  to  the  complex  modulus  of  the 
ATF  model,  equation  (35),  rewritten  here  for  convenience: 
fd  -  A)  +  i(G>/B)l 


E* (w)  =  E 


(62) 


1  +  i(w/B) 


Setting  the  right  hand  side  of  equation  (61)  equal  to 
equation  (62)  results  in  the  comparison  equations  which 
reveal  values  which  the  ATF  model  parameters  must  assume  in 
order  to  model  thermoelastic  damping. 

Comparison  of  the  denominators  of  the  two  moduli 
results  in  the  following  equation: 


B  =  (Kki )/Cv 


(63) 
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The  real  parts  of  the  moduli 's  numerators  results  in  the 
relationship 

E(1  -  4 )  =  Et  (64) 

This  equation  shows  that  the  modulus,  E,  of  the  ATF  model, 
when  used  to  describe  thermoelasticity,  is  related  by  an 
expression  in  A  to  the  isothermal  modulus,  Et . 

Comparison  of  the  imaginary  parts  of  the  numerators  of 
the  complex  moduli  resulting  from  the  thermoelastic  and  ATF 
models  gives  the  following  relation: 

E/B  =  (EtCv  +  0<T2ToET2)/(A<fk2)  (65) 

Substituting  the  expression  for  E  found  from  equation 
(64),  and  the  expression  for  B  from  equation  (63)  results  in 
(EtCv)/((1  -Zl)/^kn  =  (EtCv  +  aT^ToETM/('«^k2)  (66) 

Solving  this  for  A  results  in  the  value  of  A  in  terms  of  the 
thermoelastic  parameters: 

A  =  (EtO*t  ^To  )  /  (EtOCt  ^To  +  Cv)  (67) 

Solving  equation  (64)  for  E  with  the  value  of  A  from 
equation  (67)  gives 

E  =  Et/(1  -  A)  =  [(ETOfT)2To  +  EtCv]/Cv  (68) 

The  ATF  model,  then,  reproduces  the  thermoelastic 
model,  given  that  the  parameters  of  the  ATF  model  are  as 
fol lows : 

1/B  =  Cv/(/Vk2) 

A  =  (EtO*t  ^To  )  /  (EtO*t  2To  +  Cv) 

E  =  Et[1  +  (EtToOt^ )/Cv]  (69) 
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It  can  be  shown  (Torvik,  1990;  Appendix  B)  that  the  E 
of  equation  (69)  is  precisely  the  Young's  modulus  for  an 
adiabatic  deformation.  Thus,  E  =  E»  ,  and  from  equation 
(64) 

A  =  1  -  Ei/Eft  (70) 

Invoking  the  definition  of  A  ,  A  -  6^/(Ect)  ,  gives 

=  a(EA  -  Et)  (71) 

Therefore,  the  free  energy  defined  by  equation  (22)  is 
seen  to  require  that  the  coupling  parameter,  5,  of  equation 
(22)  be  (5  =  [a(EA  -  Et)]^/^  .  Likewise,  the  parameters  a 

and  E  of  equation  (22)  are  the  coefficient  of  thermal 
expansion  and  adiabatic  modulus,  repectively. 

Comparison  to  the  Viscoelastic  Model.  Recall  that  the 
generalized  linear  constitutive  equation  for  a  viscoelastic 
material  is  given  by  equation  (8),  repeated  here  for 
clarity : 

M  N 

^  bad^q/dt"  =  ^  andi^e/dt"  (72) 

m=0  n=0 

This  can  be  written  in  the  simplest  case  of  viscous  damping 
as 

<T(t)  =  ao6(t)  +  aid6/dt  (73) 

in  which  bo  is  set  to  unity  without  loss  of  generality.  In 

this  case,  ao  is  the  elastic  modulus  Eo .  Using  this  in 

equation  (72)  results  in 

M  N 

Yj  bad-a/dt-  +  (7(t)  =  E„€(t)  +  2Iandne/dtn 
m=l  n=l 


(74) 


Taking  the  Fourier  transform  of  equation  (74),  where 


temporal  derivatives  introduce  (iu)),  with  a{u)  and  €(u>) 
being  the  Fourier  transforms  of  a(t)  and  e(t),  respectively, 
and  dividing  both  sides  of  the  result  by  6(ci>)  results  in  an 
equation  for  the  complex  modulus: 

N 

o(<jj)  Eo  +  r  an(iw)^ 

E*((j)  =  -  =  - -  (75) 

€{u))  1  +  X  b,(ia>)» 


The  complex  modulus  for  the  Kelvin  solid  includes  only  one 
term  from  the  summation  of  the  numerator,  and  the  complex 


modulus  becomes 


E*(a})  =  Eo  +  i(aiw)  (76) 

This  model  of  viscoelasticity  increases  linearly  with 
frequency.  In  turn,  the  model  of  a  Maxwell  fluid  includes 
only  one  term  from  the  summation  in  the  denominator  of 
equation  (75),  with  the  complex  modulus  given  by 

E’(u))  =  Eo/(l  +  i(biu»)  (77) 


which,  after  multiplication  of  the  complex  modulus 
denominator,  gives 

E*(u))  =  Eo(l  -  ibiw)/(l  +  b2w2) 

Taking  one  term  from  the  summations  of  each  of  the 
and  denominator  of  equation  (75)  results  in 
Eo  +  iwai 

E*(w)  =  - 

1  iu)b  ] 


of  the 

(78) 

numerator 


(79) 


Comparison  of  equation  (79)  to  equation  (35),  the 
equation  of  the  complex  modulus  for  the  single-ATF  model, 
results  in  the  following  equations: 


(80) 


1/B  =  bi 
E(1  -  4)  =  Eo 
E/B  =  ai 

This  comparison  is  easily  made,  and  the  ATF  model 
reduces  to  the  viscoelastic  model  with  3  parameters  —  Eo , 
ai ,  and  b] ,  if  the  constants  are  equated  as  follows: 

B  =  1/bi  (81) 

E  =  ai/b,  (82) 

d  =  1  -  (Eobi)/ai  (83) 

Inspection  of  equation  (79)  shows  that  the  low 
frequency  modulus  must  correspond  to  the  isothermal  elastic 
(rubbery)  modulus.  En  is  identified  as  such.  Conversely, 
as  the  frequency  goes  to  infinity,  the  real  part  of  the 
modulus  becomes  ai/bi.  Thus,  this  ratio  is  the  adiabatic 
(glassy)  modulus.  Accordingly,  the  E  used  in  the  definition 
of  the  free  energy  must  be  the  adiabatic  modulus  and  the 
parameter  4  is 

4=1-  Er/Ei  (84) 

Recall  from  equation  (36)  that  4  =  ^^/(Ea)  ,  so  the 

coupling  parameter,  <5,  must  be 

<5  =  [a(EA  -  Et)]  1/2  (85) 
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III.  Fractional  Derivatives  In  the  Description 
of  Material  Damping 

Capitalizing  on  the  previously  made  observation  that 
the  use  of  fractional  derivatives  in  material  models  enables 
the  description  of  material  damping  over  a  much  broader 
frequency  range,  this  research  effort  proposes  to  improve 
upon  the  accuracy  of  Lesieuture's  model  over  a  broader 
frequency  range  through  the  use  of  fractional  derivatives. 
Coupled  material  constitutive  relations  will  be  developed 
using  the  concept  of  augmenting  thermodynamic  fields,  with 
non-integer  differentials  allowed  in  the  resulting  partial 
differential  equations. 

Viscoelasticity  Using  the  4-Parameter  Model.  Bag ley  and 
Torvik  have  found  the  use  of  stress-strain  laws, 
constitutive  relationships,  or  equations  of  state  which 
employ  fractional  or  generalized  derivatives  to  be  valid  as 
a  material  damping  model  (Torvik  and  Bagley,  1987). 
Relationships  using  a  generalized  derivative  have  been  shown 
to  have  a  sound  theoretical  basis,  and  are  effective 
descriptors  of  the  dynamic  behavior  of  real  materials.  A 
brief  review  of  Torvik  and  Bagley' s  work  will  be  presented 
here . 

Fractional  Calculus.  Ross  presents  a  basic  history  of 
the  fractional  calculus,  in  which  derivatives  and  integrals 
of  fractional  order  are  defined,  and  states  that  it  is 
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nearly  as  old  as  the  calculus  of  integer  orders  (Ross, 

1977).  The  extended  Riemann-Liouvi 1 le  definition  of  the 
generalized  derivative  is  given  as  (Oldham  and  Spannier, 
1974) 

dqx(t)  1  d  ^x(t-T) 

D’(x(t))  =  -  - - /  -  dr  (86) 

dri  r(l-q)  dt  o  to 

for  0  1  q  <  1  .  r  is  the  gamma  function,  defined  by 

CD 

y^(l-q)  =  e-xx-qdx  (87) 

O 

At  first  glance,  this  would  appear  to  be  a  formidable 
equation,  yet  when  one  considers  the  Fourier  (or  Laplace) 
domain,  the  generalized  derivative  manifests  itself  as  a 
fractional  power  of  u>  (or  s)  (Caputo,  1976).  The  Fourier 
transform  is  defined  as 

CO 

FT(x(t)]  =  y'x(t)e-iwtdt  (88) 

-CD 

If  x(t)  *  0  for  t  <  0  ,  then  the  Fourier  transform  can 

be  written  as 

CO 

FT[x(t)]  =  y'x(t)e-iwtdt  (89) 

o 

Taking  the  Fourier  transform  of  the  generalized  derivative 
of  equation  (86)  gives 

FT[Dfl[x(t)]]  =  (icj)qFT(x(t)]  (90) 

The  restriction  placed  on  q  is  that  it  be  a  nonnegative 
real  number  less  than  one.  For  engineering  applications,  it 
is  assumed  that  an  irrational  number  can  always  be 
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approximated  by  a  rational  number,  allowing  for  the 
restriction  of  q  to  being  a  rational  number.  Thus  only 
fractional  derivatives  of  rational  order  will  be  considered. 

Fractional  Derivatives  in  Viscoelastic  Theory.  The 
generalized  differential  operator  form  of  the  stress-strain 
constituitve  equation  is  given  by  Christensen  (Christensen, 
19B2:14)  as 

M  d»a(t)  N  dne(t) 

E  b. -  +  <T(t)  =  E€(t)  +  E  an -  (91) 

m=l  dt«  n=l  dtf' 

The  Fourier  transform  of  equation  (91)  yields 

W  _  _  N  _ 

E  bn  (icj)  »^(u>)+a(aj)  =  Ei(u))+ 53  an  ( iu))  r'f  (u»)  (92) 

m=l  n=l 

Through  manipulation  of  this  equation,  a  complex  modulus  is 

formed : 

N 

a(u))  E  +  ar,  (iu))'^ 

E*(u»  =  I -  =  - fi -  (93) 

f(uJ)  1  E  b,(iu>)* 

Now,  if  fractional  derivatives  are  allowed  in  the 
general  form  of  the  relationship  for  linear  viscoelasticity, 
the  result  is  (Bagley  and  Torvik,  1979) 


H  N 

E  b.Dq  [(7(t)]+<T(t)  »  Eo€(t)+E  EnDq[e(t)]  (94) 

m=l  n=l 


where  q  is  a  number  between  zero  and  one. 

This  is  the  form  appropriate  for  a  one-dimensional 
description  of  axial  deformation.  Many  materials  can  be 
accurately  modeled  by  replacing  each  sum  in  equation  (94)  by 


a  single  term  involving  a  fractional  derivative 
<T(t)  +  bD<J[(7(t)]  =  Eo€(t)  +  EiDq[€(t)] 


(95) 


where  b,  q,  Eo,  and  Ei  are  the  real  parameters  of  this  four 
parameter  model  (Bagley  and  Torvik,  1986).  The  Fourier 
transform  of  this  model  yields  a  complex  modulus  of  the  form 
Eo  +  E\(iciJ)<3 

E*(u})  =  -  (96) 

1  +  b(ia})r 

with  Eo  and  Ei/b  being  the  rubbery  and  glassy  moduli, 
repectively . 

Thermodynamic  constraints  require  the  following 
conditions  of  the  parameters  of  equation  (96): 

Eo  >  0  El  >  0 

b  >  0  Ei/b  >  0  (97) 

These  constraints  insure  a  themodynamical ly  well-behaved 
fractional  calculus  model  of  viscoelastic  behavior  (Bagley 
and  Torvik,  1986). 

Damping  Ratio.  In  order  to  solve  for  the  damping 
ratio,  equation  (96)  must  be  further  simplified  using 
equations  for  complex  numbers  found  in  Appendix  A.  Using 
these  simplifications,  equation  (96)  may  be  rewritten  as 

Eo  +  El  (uj) “J  [cos(q7r/2)  +  i  8in(q7r/2)] 

E*(a»  =  -  (98) 

1  +  b(w)  1  (cos (q7r/2)  +  i  sin(q77‘/2)] 

Letting  x  =  q7r/2  and  multiplying  equation  (98)  by  the 
complex  conjugate  of  the  denominator  results  in  the 
following  for  the  complex  modulus: 

E'{w)  -  N/D  (99) 

where 
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N  =  Eo  +  (Eob  El  )a><icos(x)  +  EibtjZ'J 


+  i  [(El  -  Eob)u>qsin(x) ]  (100) 

and 

D  =  1  +  2ba;<Jcos(x)  +  (101) 

Recall  from  equation  (4)  that  the  loss  factor,  t]  , 
equals  the  Imaginary  part  of  the  complex  modulus  divided  by 
its  real  part: 

(El  -  Eob)tj<Jsin(x) 

77  =  -  (102) 

Eo  +  (Eob  +  El  )tJ‘>cos(x)  +  Eibu>2<5 

with  the  damping  ratio  being  one-half  the  loss  factor. 

A  plot  of  the  loss  factor  versus  the  frequency  of 

equation  (102)  is  included  as  figure  4.  Note  the  strong 

frequency  dependence  of  the  damping  ratio.  Also, 


Frequency 

Figure  4.  Loss  Factor  vs  Frequency  of  4-Parameter 
Fractional  Derivative  Model 


33 


examination  of  equation  (97)  shows  that  the  thermodynamic 
constraints  insure  that  the  loss  factor  is  positive  for  all 
frequencies  (Bagley  and  Torvik,  1986).  Values  used  for  the 
plot  of  figure  4  are 

Eo  =  5.04  X  lOb  N/m2 

El  =  2.27  X  105  (N  s)/m2 
b  =  2.80  X  10-4  s 
q  =  0.64 

These  values  correspond  to  Nitrile  Rubber  1479  resulting 
from  a  "least-squares"  fit  to  data  previously  obtained 
(Bagley  and  Torvik,  1986:137). 

Single-ATF  Model  With  Frational  Derivatives.  Torvik  and 
Bagley  (Torvik  and  Bagley,  1987:125)  noted  that  most 
materials  used  in  damping  application  are  polymers  showing  a 
strong  frequency  dependence.  The  general  viscoelastic 
relationships  required  to  describe  and  adequately  model  this 
dependence  usually  require  a  large  number  of  terms.  They 
have  noted  that  the  frequency  dependence  of  real  materials 
is  more  closely  approximated  by  models  with  fewer  terms,  but 
whose  material  constitutive  time  derivatives  are  of  a 
fractional  order.  With  four  and  five  parameter  models,  they 
have  accurately  modeled  material  damping  behavior  over  broad 
ranges  of  frequency,  some  even  approaching  eight  decades. 

Although  Lesieutre  declares  this  approach  as 
"cumbersome,"  that  is  not  the  case  and  the  ease  of  its  use 
in  application  to  the  ATF  model  will  be  shown  here. 
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Constitutive  Equations.  The  equation  for  the  Helmholtz 
free  energy  is  as  was  used  before,  and  can  be  found 
as  equation  (22).  The  fractional  derivative  will  be  applied 
to  the  equilibrium  equation  of  irreversible  thermodynamics 
as  follows  from  equation  (26): 

a<J^/at<i  -  -B(^  -  1)  =  +  hd€/a  (103) 

Here,  we  assume  the  rate  process  governing  the  restoration 
of  the  equilibrium  state  is  a  process  of  fractional  order  in 
time.  This  equation  together  with  the  balance  of  momentum 
relationship,  equation  (30),  gives  two  coupled  partial 
differential  equations: 

pii  -  Eu"  =  0 

+  B4  -  Biu'/a  *  0  (104) 

In  the  Fourier  domain  these  equations  become  functions 
of  u).  Taking  the  Fourier  transform  of  equation  (104)  with 
time  derivatives  introducing  iu)  and  spatial  derivatives 
introducing  ik  gives 

P(iw)?  U  -E(ik)2  U  +  d(ik)  Z  «  0  (105) 

(iu>)<)  Z  +  B  Z  -  B(5(ik)/a  U  *  0  (106) 

Solving  equation  (106)  for  Z  results  in 

Z  ■  IBk^/raB  +  a(iw)<))  u  (107) 

Substituting  equation  (107)  into  equation  (105)  and 
canceling  the  U  gives  an  equation  that  can  be  algebraically 
manipulated  to  give  the  complex  modulus: 

Ek2  -  Bk2()2/[otB  •♦■a(ia>)'>]  «  0  (100) 
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Recall  (Appendix  A)  that  {iu^)'3  can  be  written  as 
(io^lq  ^co<i  [cos(q7r/2)  +  i  8in(q7r/2)]  (109) 

Substituting  x  =  q7r/2  into  equation  (109)  and  the  result 
into  equation  (106)  results  in 

-pu)2  +  Ek2  -  2 /{a(B+u;<Jcos(x)-iw‘>sin(x) ]  >  =  0  (110) 

Algebraic  manipulation  of  equation  (110)  results  in  a 
familiar  form; 


k2 


1  - 


1  +  (u><J/B)  [cos(x)  +  i  sin(x)] 


(111) 


where 

A  =  (5V(Ea) 


(112) 

The  right  hand  side  of  equation  (111)  is  recognized  as  the 
complex  modulus,  which  can  be  rewritten  as 

1  -d+  (uJ'i/B)  (cos(x)  +  1  sin(x)] 


E-(u))  =  E 


(113) 


1  +  (u>=3/B)  [cos(x)  +  i  sin(x)) 

Multiplying  the  numerator  and  denominator  of  the  complex 
modulus  of  equation  (113)  by  the  complex  conjugate  of  its 
denominator,  1  +  (wj /B) cos (x)  -  i(u>‘J/B)sin(x)  ,  results  in 
the  following  after  some  algebraic  manipulation: 

r 

1  (cj'J /B)cos(x)  -  i  (a>'>/B)sin(x)| 


E*  (u) 


1  -A 


(114) 


1  +  2(u;i/B)cos(x)  + 

Damping  Ratio.  The  loss  factor  can  be  found  again  from 
equation  (4): 

A  (cj*>/B)sin(x) 


l  +  2((p'l/B)C08(X)  +  (u;<I/B)  2-dU'*'(ci>’/B)co8(x)l 


(115) 


Two  assumptions  are  again  made  to  arrive  at  the  damping 
ratio.  They  are:  (1)  small  damping,  which  in  turns  implies 


that  the  deunping  ratio  is  half  the  loss  factor,  and  (2) 
d  <<  1  .  which  will  cancel  all  terms  involving  d  from  the 

denominator.  This  results  in  an  equation  for  the  damping 
ratio  in  terms  of  the  frequency: 

((5V2Eof)  (a><J/B)  sin(q7r/2) 

4  - - (116) 

1  +  2(tj<i/B)  cos(q7r/2)  +  (fj'J/B)^ 

By  setting  q  of  equation  (115)  equal  to  1  the  damping 
ratio  reduces  to  that  of  Lesieutre's  single-ATF  model,  since 
cos(q7r/2)  goes  to  zero  and  sin(q7r/2)  goes  to  one.  The 
fractional  derivative  single-ATF  model  can  then  be 
considered  a  more  general  case  of  Lesieutre's  model. 

Figure  5  is  a  plot  of  the  damping  ratio  versus 
frequency  for  three  values  of  q  in  the  fractional  derivative 


Figure  5.  Single  ATF  Model  with  Fractional  Derivatives 


single-ATF  model.  Values  for  the  other  parameters  were  as 
fol lows : 


B  =  1.0 

d  =  0.1 


As  the  value  of  q  decreases  the  curve  flattens  out,  as  well 
as  opening  up. 

Complex  Modulus  Formulation.  An  alternate  method  of 
arriving  at  the  loss  factor,  and  thus,  the  damping  ratio, 
might  be  termed  the  "complex  modulus"  formulation.  Starting 
with  equations  (24)  and  (103)  results  in,  after  the  Fourier 
transform: 


d(oj)  =  -  d2(u)  (117) 

=  -B2(w)  +  ih6/a)€(cj)  (118) 

where  5(u>)  is  the  Fourier  transform  of  a(x,t),  e(u))  is  the 
Fourier  transform  of  €(x,t),  and  Z(w)  is  the  Fourier 
transform  of  ^(x,t).  Solve  equation  (118)  for  Z(u)): 

Z(w)  =  {B  (5  /  (OB  +  a(iw)<i]>  €  (a;)  (119) 

Substituting  the  result  into  equation  (117)  and  dividing  by 
€(u?)  results  in  an  equation  for  the  complex  modulus: 


a(cj)/€{uj)  *  E  -  B()2/[o‘B  +a(i<*;)'i] 

Recalling  ii  «  cos(q7r/2)  i  sin(q7r/2)  and  letting 
X  B  q7r/2  ,  equation  (120)  becomes 


(120) 


aiu) 

€  (<J) 


1  -A 


l  +  (u'*/B)cos(x)  +  i  (iL*<i/B)sin(x) 


(121) 


where  d  is  defined  as  before  in  equation  (112).  This  is 
the  complex  modulus.  After  multiplying  the  numerator  and 
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denominator  of  the  Innermost  term  by  the  denominator's 
complex  conjugate,  equation  (121)  takes  on  the  familiar  form 
of  equation  (114): 


l  +  (tj'>/B)cos(x)-i  (a><J/B)sin(x 


E'(to) 


E 


1  -  A 


(122) 


[  1^1  +  2(tj9 /B)cos(x)  +  (tj<J/B)2j| 

Derivation  of  the  damping  ratio  from  this  point  follows  the 
previous  derivation. 


Comparison  of  the  ATF  model  with  Fractional  Derivatives  to 
the  4-Parameter  Model.  Recall  from  equation  (98)  that  the 
'expression  for  the  complex  modulus  of  the  four  parameter 
model  of  viscoelasticity  is 

Eo  +  Eia)r(cos(y)  +  i  sin(y)] 

E*(u))  =  -  (123) 

1  +  bu;''(cos(y)  +  i  sin(y)] 

where  y  =  rn/2  from  the  derivation  of  equation  (98). 

Equation  (123)  can  now  be  compared  to  equation  (113), 
the  equation  of  the  complex  modulus  for  the  single-ATF  model 
with  fractional  derivatives.  The  following  equations  of 
comparison  result: 

1  +  (uJ'J/B)  (cos(x)  +  1  sin(x)] 


=  1  +  (bu)*^ )  (cos(y )  +  i  sin(y)) 

(124) 

E(1  “  d )  »  Eo 

(125) 

E(u>i/B)  [cos(x) 

+  i  sin(x)] 

■  (Eiw )  [cos(y)  +  i  sin(y)] 

(126) 

Examination  of  all  equations  show  that  r  must  equal  q 
in  order  for  the  terms  involving  w  and  the  trigonometric 


functions  to  equate.  Therefore,  the  first  equivalency 
equation  is 

r  =  q  (127) 

which  in  turn  implies  x  =  y 

Substituting  equation  (127)  into  equation  (124)  and 
solving  for  B  shows  that  the  second  condition  for  model 
equivalency  is 

B  =  1/b  (128) 

Using  equation  (128)  in  equation  (126)  gives  the 
relationship: 

E/B  =  El  (129) 

The  final  relationship  of  comparison  is  that  of 
equation  (125),  repeated  here  for  clarity: 

E(1  -  4)  =  Eo  (130) 

Equations  (127)  through  (130)  are  the  equations  that 
show  equivalency  between  the  4-Parameter  model  and  the  ATF 
model  with  fractional  derivatives.  It  should  be  noted  that 
no  assumptions  were  made  as  to  small  4.  Recall  that  the 
glassy  modulus  equals 

Eo  =  Ei/b  =  (E/B)/(l/B)  =  E  (131) 

Comparing  this  to  the  equation  for  the  rubbery  modulus, 
equation  (131)  shows  that  the  ratio  of  the  glassy  modulus  to 
the  rubbery  modulus  is 

Eq/Eo  =  E/(E(1  -  4)]  (132) 

In  many  materials  the  glassy  modulus  is  three  orders  of 


magnitude  larger  than  the  rxibbery  modulus.  This  being  the 
case,  A  must  equal  0.999  or  greater. 

As  can  be  seen,  A  is  not  small  as  compared  to  one  in 
this  case.  Bagley  and  Torvik  have  used  the  4-Parameter 
model  to  produce  curves  with  accurate  results  over  a  very 
broad  frequency  range  for  many  materials,  principally 
polymers.  The  assumptions  that  d  <<  1,  then,  is 
inapproriate  for  use  with  4-parameter  viscoelastic  modeling. 

For  comparison  purposes,  the  plot  of  the  loss  factor 
versus  frequency  for  the  single-ATF  model  without  the  small 
^  assumption  is  shown  as  figure  6.  The  values  used  for  the 
plot  are  as  follows: 

B  =  1 

d  =  0.9995 


Figure  6.  Single-ATF  Model  Without  Small  d  Assumption 


No  cancellation  of  A  from  the  denominator  of  the  loss  factor 
was  made  in  producing  figure  6.  Figure  7,  for  comparison 
purposes,  is  a  plot  of  the  loss  factor  using  the  same  values 
of  B  and  but  with  denominator  terms  involving  d  canceled. 


Figure  7.  Single-ATF  Model  With  Small  d  Assumption 


IV.  Two-ATF  Model 


Lesieutre's  development  of  the  model  for  the 
one-dimensional  vibration  of  an  isotropic  rod  was  selected 
for  ease  of  illustration  of  the  ATF  modeling  method.  He 
states  that  in  practice  additional  augmenting  thermodynamic 
fields  can  be  used  as  needed  to  better  approximate 
experimental  data  over  the  frequency  range  of  interest 
(Lesieutre,  1989:19).  In  order  to  illustrate  the  ease  or 
difficulty  of  adding  additonal  ATFs ,  the  following 
developments  of  two-ATF  models  were  performed. 

Two-ATF  Model  with  Integer-order  Derivatives.  For  the 
development  of  the  two-ATF  model  with  integer-order 
derivatives  a  second  augmenting  thermodynamic  field  is  added 
to  the  free  energy  equation.  The  formulation  of  the  damping 
ratio  followed  that  of  Lesieutre's  single-ATF  model,  with 
the  addition  of  equations  involving  the  second  ATF. 

Constitutive  Equations.  Once  again,  the  Helmholtz  free 
energy  is  employed  to  arrive  at  the  const ituitive  equations. 
We  hypothesize 

f  =  1/2E6Z  +  l/2ai^i2  + 

-  ~  ~  (133) 

Here,  £  is  the  strain  field,  and  are  the  two  augmenting 
thermodynamic  fields,  <5i  and  d?  are  the  strengths  of  the 
coupling  between  strain  and  and  ^2,  respectively,  and 
is  the  strength  of  the  coupling  between  the  two  ATFs,  ^1  and 


The  stress  field  is  again  found  from  the  relation 
o  =  df/^e  =  E€  -  -  67^?  (134) 

and  the  two  affinities  A\  and  A?  can  be  found  using  equation 
(25),  which  results  in 


Ai  =  -at/a^i  = 

+0*3^2 

(135) 

A?  =  ~af/a^2  =  67^ 

-  ct7  ^7  + 

(136) 

Assuming  the  irreversible  thermodynamic  processes  are 
described  by  first  order  rate  equations,  the  time  rate  of 
change  of  ^i,  or  is  assumed  to  be  proportional  to  the 

deviation  from  an  equilibrium  value: 

a^/at  =  -Bi(^;  -  ^i)  i  =  l,2  (137) 

with  the  value  of  found  from  the  value  of  ^ ,  at 
equilibrium,  or  when  the  affinity  is  zero.  Here  it  is 
assumed  that  there  is  no  coupling  between  the  two  ATFs  in 
the  equilibrium  relationship  of  equation  (137),  which  would 
be  indicated  by  an  equation  of  the  form 

a^/at  =  -  |i)  -  Ci(^2  -  (2)  (130) 

for  the  first  ATF,  with  a  similar  equation  for  the  second 
ATF.  Such  interaction  are  not  used  in  the  developments  of 
the  two-ATF  models. 

Setting  Ai  in  equation  (135)  equal  to  zero  and  solving 
for  4 1  gives 

=  6i€/0l\  +  Ct3^2/fl»i  (139) 

Similarly,  for  ^2. 

^2  =  67^/017  +  Ol3^i/Ol2  (140) 

Substituting  these  results  into  equation  (137)  gives,  for 


the  time  dependence  of  the  ATFs,  two  equations: 

a^i/at  =  -  (5i€/ai  -  a3^?/aj)  (141) 

d^^/dt  =  -  bii/ci'i  ~  (142) 

Notice  that  now  and  are  coupled  in  and  ^7.  as  well 
as  in  strain,  €  . 

Balance  of  linear  momentiun  requires  that,  for  a 
one-dimensional  deformation, 

da/dx  -  pa2u/at2  =  0  (143) 

The  derivative  of  equation  (134)  with  respect  to  x  gives 

acr/ax  =  Eae/ax  -  6xdU/dx  -  6 76^7 /ex  (141) 

•which,  when  substituted  into  equation  (143)  and  recalling 
€(x,t)  =au(x,t)/ax  =  u'  ,  results  in  a  partial 
differential  equation  coupled  in  u,  and 

pu  -  Eu"  +  +  67U'  -  0  (145) 

where  the  dots  represent  time  derivatives  and  the  primes 
represent  spatial  derivatives.  This  equation  together  with 
equations  (141)  and  (142)  are  the  governing  equations  for 
longitudinal  oscillations  of  an  axial  rod. 

Damping  Ratio.  In  order  to  solve  the  constitutive 
equations  with  greater  ease,  the  Fourier  transform  is 
Invoked,  with  time  derivatives  Introducing  iu),  and  spatial 
derivatives  Introducing  ik,  which  results  in  the  following 
equations : 

p(iw)2U  -  E(ik)2U  +  di(ik)Zi  +  a2(ik)Z?  =  0 
Bi($i(ik)U  —  011(85  iu))Zi  +  BiOtsZ?  =  0 
B2<5?(ik)U  +  BrOsZi  -a2(B?  +  iu))Z?  =  0  (146) 
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These  equations  may  be  expressed  in  matrix  form  as 

-pu>2  +Ek2  i^ik  idik  U 

iB]d]k  “0*1  (B]  +  icj)  B]0(3  '  Zi  ’  =  {0} 

iB2(52k  B20*3  -02(82  +  iw)  Z2  (147) 

' .  ) 

This  is  a  linear  system  of  equations  whose  solution  is 
determined  by  finding  the  determinant  of  the  matrix  and 
setting  it  equal  to  zero: 

(“pw2  +  Ek2 )  [  (-0*1 )  (Bi  +  iw)  (-0*2 )  (B2  +  iw)  -  BiB20*32] 

-  ik(5i  [iBik(5i  (-0*2 )  (B2  +  iw)  -  iBi B2ko*3(5  2  ] 

+  ik(52  [iBik^  1  (B20<3  )  +  iB2ko*ia2(Bi  +  iw)]  =  0  (148) 

Carrying  out  the  algebra  and  rearranging  equation  (148) 

gives 

(-pu,2  +  Ek2 )  [r  -  wV(BiB2)  +  iw/Bs] 

-  k?E[Ji  +  ^2  +  43  +  iw(4i/B2  +  42/Bi)]  =  0  (149) 

where 

y  =  1  -  032/(0102) 

di  =  (5i2/(Eai) 

d2  =  (J22/(E02) 

d3  =  2  ((5  1  (5  203  )/(Ea, 0*2) 

l/Bs  =1/Bi  +  I/B2  (150) 

Dividing  equation  (149)  by  7  -w2/{BiB2)  +  iw/Bs  , 
dividing  the  result  by  p,  and  factoring  k^/p  from  this 
results  in  an  equation  from  which  the  complex  modulus  can  be 
derived : 

k2  E[4i  42  43  +  ia>(4i/B2  +  42/Bi)] 

-<j2  +  —  E - =  0  (151) 

P  y  -  U2Z/(BiB2)  +  iCu>/B3 
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The  bracketed  cfuantity  of  equation  (151)  is  the  complex 
modulus  E*(w),  and  the  loss  factor  r) .  therefore,  can  be 
found  from  equation  (4) .  Multiplying  the  complex  modulus  by 
the  complex  conjugate,  y  -u)2/(BiB2)  -  iw/Ba  ,  of  the 
denominator  of  the  second  term  within  the  brackets  of 
equation  (151)  results  in  the  following  expression  for  the 
complex  modulus: 


E*  (cj) 


E 


Nre  +  i  Ni. 

(>■  -  tJV(BlB2)]2  +  (CJ/B3)2 


(152) 


where 


N 


r  e 


and 


[>’~4j2/(BiB2)][>'-tt/2/(BiB2)  +^i  +^2  +4’] 

+  (<*>/B3  )  (4it*^/B2  +  A7CJ/B1  +  tJ/Ba)  (153) 


Ni.  =  [y  -  (jV(BiB2)]  [4ia;/B2  *  A7CJ/B1] 

-  ((J/Ba)(4i  +4?  +  4a)  (154) 

The  loss  factor,  found  using  equation  (4),  is  given  by 


r]  =  Im{E’ (a2)}/J?e{E*  (<*>)} 

=  Ni./Nre  (155) 

where  Nj.  and  Nre  are  given  by  equations  (154)  and  (153), 
respectively. 

With  the  assumptions  that  damping  is  small,  and  that 
<<  1,  A  7  <<  1,  and  da  <<  1,  the  damping  ratio  then 


becomes 


f- 


(tj/Ba )  (4i +42  +  43 )“  /(BjB? )  1  (a»4i  /B2+<*>42  /Bi  ] 
2{(y  -  <J2/(BiB2)]  2  +  (W/B3)2> 


(156) 
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In  order  to  check  that  this  reduces  to  the  case  when 


only  one  ATF  was  used,  set  as  =  0  and  d?  =  0,  which  from 
equation  (150)  assures  A?  =  ds  =  0  and  y=  1.  Equations 
(153)  and  (154)  then  become 

Nrel  =  1  -  (2a>2)/(BiB2)  +a>V(BiB2)2 

♦  ((jMi)/(B2B3)  +W2/B32  (157) 

Ni.i  =  Aicj[1/B2  -  -  I/B3]  (158) 

and  the  loss  factor  is  an  equation  in  terms  of  di ,  Bi ,  and 
B2,  with  no  coupling  between  the  second  ATF,  ^2,  and  the 
strain,  €,  or  the  first  ATF,  ^1.  The  loss  factor  is  then 
^  =  Ni.i/Nrei 


6J/I1  (Bi+B2)/(BiB2)  -  [1  -  w  V(BiB2))  [tJdi/B2] 

=  -  (159) 

(1-w2/(BiB2)1  2  +  (6>(Bi+B2)/(BiB2)]  2 

As  can  be  seen,  equation  (159)  is  not  equivalent  to  twice 

equation  (41),  the  loss  factor  for  a  single  ATF,  in  that  the 

material  constant  B2  is  still  present.  Taking  the  limit  of 

equation  (159)  as  B2  goes  to  infinity  results  in 


cjdi  /Bi 

r]  =  -  (160) 

1  +  (j  ^  /B 1  2 

Examination  of  equation  (160)  shows  it  to  be  identical 
to  the  value  of  the  loss  factor  for  the  case  where  only  one 
ATF  was  employed.  The  loss  factor  of  the  two  ATF  model, 
therefore,  reduces  to  that  of  the  one  ATF  model  exactly. 

Characteristics  of  the  Two-ATF  Model.  Figure  8  is  a 
plot  of  the  loss  factor  versus  frequency  for  the  two  ATF 
model  with  the  values  of  the  parameters  as  follows: 


48 


Bi  =  1 
B2  =  10 

Ai  =  0.1 

A?  =  0.05 


Frequency 

Figure  8.  Loss  Factor  vs.  Frequency  -  2-ATF  Model 


Here  Bi  and  B?  have  the  dimensions  of  frequency,  and  Ai  and 
A?  are  dimensionless.  These  values  were  selected  only  to 
allow  for  ease  in  comparison  as  the  parameters  of  the  model 
are  changed.  Again,  0(i  and  are  assigned  values 
equivalent  to  the  magnitudes  of  Bi  and  B?,  respectively 
(Lesieutre,  1989:48,  146). 

It  will  be  noted  that  As  is  a  function  of  <5i  and  67.  as 
well  as  «! ,  ot2 ,  and  OI3 .  From  equation  (150) 

(5i  =  (AiEOfi)i/2  (161) 

67  =  (A2Ea2)W2  (162) 
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4  q  = 


Eotioc: 


OflO? 


oeioe? 


(163) 


a3 ,  then,  can  be  found  from  a  relationship  of  ai  and  <32: 


03  =  C’  (otia? )  1  /  2 


(164) 


43  =  C*  (4i42)i/2 


(165) 


where  C*  is  a  constant  parameter  that  can  be  varied  as  the 
other  material  parameters.  For  the  plot  in  Figure  8  C’  was 
set  equal  to  0.1,  and  ds  =  X  10- 2 

Figure  8  shows  the  "double  hump"  characteristic  of  the 
loss  factor  when  two  ATFs  are  employed,  with  the  peak  of  the 
first  maximum  at  a  frequency  equal  to  the  value  of  Bi  and 
the  second  maximum  falling  on  a  frequency  equal  to  the  value 
of  B? .  Notice  that  two  maxima  are  discernable,  although 
they  almost  merge  into  a  single  curve.  This  discernabi 1 ity 
is  not  clearly  the  case  in  some  of  the  later  plots. 

Figures  9  through  18  are  plots  of  the  loss  factor 
versus  frequency  with  different  parameters  of  the  two-ATF 
model  varied  for  each  plot. 

Figure  6  shows  the  loss  factor  versus  frequency  with 
B]=0.2.  The  first  maximum  has  shifted  to  the  left  and 
corresponds  to  a  frequency  of  0.2.  The  two  maxima  can  be 
seen  as  separate  maxima.  This  is  also  obvious  in  figure  10, 
in  which  the  value  of  B?  is  increased  to  90.  Again,  the 
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Frequency 


Figure  9.  Two-ATF  Model  with  Bi=0.2 
(B?=10,  Ai=0.1.  A?=0.05,  A3=v^  X  10-2) 


Frequency 


Figure  10.  Two-ATF  Model  with  B2=90 
(Bi=l,  Aj=0.1,  A?=0.05,  A3='/r  X  10-2) 
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second  maximum  shifts  to  a  value  corresponding  to  90. 
Decreasing  the  value  of  Bi  also  tends  to  decrease  the 
maximum  value  of  the  loss  factors  for  both  peaks  of  the 
curve,  which  also  occurs  with  an  increase  in  B2 . 

The  opposite  effect  is  apparent  as  Bi  is  increased  to 
9.0  and  B?  is  decreased  to  2  in  figures  11  and  12, 
respectively.  As  the  values  of  Bi  and  B2  approach  each 
other,  the  double  hump  collapses  into  one  maximum,  and  the 
value  of  the  loss  factor  tends  to  increase.  In  both  cases 
the  value  of  the  maximum  is  the  same,  with  the  curve 
appearing  to  be  the  same,  only  shifted  along  the  x-axis. 


Frequency 


Figure  11.  Two-ATF  Model  with  B  =9.0 
(B2  =  10,  Ai=0.1.  ^2=0.05,  d3»'/5  X  10-2) 
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Frequency 


Figure  12.  Two-ATF  Model  with  B2=2.0 
(Bi=l,  Ai=0.1,  A2=0.05,  A3=v^  X  10-2) 


Figures  13  and  14  show  the  effects  of  varying  <  with 
figure  13  corresponding  to  an  increase  cf  to  0.9,  and 
figure  14  to  a  decrease  to  0.02.  An  increase  of  tends  to 
collapse  the  separate  maxima  into  one  peak  and  to  increase 
the  value  of  the  loss  factor  at  that  maximum.  The  value  of 
the  maximum  loss  factor  is  almost  an  order  of  magnitude 


higher.  Although  barely  perceptible,  an  inflection  point, 
or  local  maximum,  occurs  at  the  frequency  corresponding  to 
B: .  With  a  decrease  in  Ai,  however,  a  decrease  in  the 


maximum  value  corresponding  to  Bi  occurs.  This  decrease  is 
not  an  order  of  magnitude  as  was  the  case  with  the  increase. 
The  value  of  the  maximum  at  u)  =  B2  decreases  as  well,  but 
remains  the  dominant  peak. 
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Figure  13.  Two-ATF  Model  with  zii=0.9 
(Bi=l,  B?=10,  d7=0.05.  d3=V5  X  10-2) 


Figure  14.  Two-ATF  Model  with  di=0.02 
(Bi=l,  B7  =  10,  d?«0.05,  d3=>/5  X  10-2) 
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Figures  15  and  16  show  similar  results  as  A7  is 
decreased  to  0.005,  with  the  peaks  merging  into  one  dominant 
maximum,  and  as  d?  is  increased  to  0.5.  In  figure  15,  the 
value  of  the  loss  factor  at  this  maximum  is  only  slightly 
less  than  that  of  figure  8,  while  the  value  of  the  loss 


Figure  15.  Two-ATF  Model  with  d? =0.005 
(Bi=l,  B2  =  10.  di»0.1,  X  10-2) 


factor  at  the  second  maximum,  barely  dlscernable  at  a 
frequency  of  10  more  as  a  deflection  point,  is  considerably 
lower.  Figure  16  shows  the  loss  factor  being  greater  over 
the  entire  frequency  range,  with  the  maximum  at  a  frequency 
of  10  being  dominant.  As  the  difference  between  di  and  d? 
increases,  the  double  maxima  merge  to  one  maximum  whose  peak 
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Figure  16.  Two-ATF  Model  with  Aj*0.5 
(Bi=l,  B7=10,  di=0.l.  X  10-2) 

is  at  the  frequency  corresponding  to  the  value  of  Bi 
associated  with  the  largest  value,  . 

As  the  strength  of  the  coupling  between  the  two  ATFs  is 
increased,  the  value  of  the  loss  factor  at  lower  frequencies 
is  increased.  Figure  17  shows  the  plot  of  loss  factor 
versus  frequency  for  da  =  \/40 ! 5  x  10- 2  ,  resulting  from 
C’'«0.9  Again,  a  maximum  at  the  value  of  frequency 
corresponding  to  B2  is  observed,  but  the  lower  maximum  has 
shifted  to  the  right  and  increased  by  an  order  of  magnitude. 
A  decrease  in  the  value  of  da  does  not,  however,  show  the 
same  trend,  as  is  evidenced  in  figure  16.  With  an  order  of 


r  requency 


Figure  17.  Two-ATF  Model  with  d,=v^40.5  x  10- z 
(B,=l,  Br=10.  di=0.1,  d?=0.05) 


Frequency 


Figure  18.  Two-ATF  Model  with  d3  =  v^0 . 5  X  10-3 
fBi=l,  B?=10,  4, =0.1,  47=0.05) 
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magnitude  decrease  in  the  value  of  C  of  equation  (164)  only 
a  very  slight  decrease  in  the  loss  factor  is  noted. 

Figure  19  shows  a  summary  of  the  effects  of  varying  the 
parameters  of  the  two-ATF  model.  The  arrow  directions 
corresponds  to  the  direction  the  curve  will  shift  with  an 
increase  of  the  parameter  shown. 


Figure  19  -  Effect  of  Parameter  Variation  on  Two-ATF  Model 

Two-ATF  Model  With  Fractional  Derivatives.  A  second 
augmenting  thermodynamic  field  is  now  added  to  the  model, 
and  derivatives  of  fractional  order  are  allowed  in  the  time 
derivative  terms  of  the  thermodynamic  equilibrium  equations. 
The  formulation  of  the  problem  follows  the  same  path  as 
previously;  however,  trigonometric  relations  are  required 
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throughout  the  algebraic  reduction  in  order  to  simplify  the 
resulting  equations  for  the  complex  modulus. 


Constitutive  Equations.  The  equation  for  the  Helmholtz 
free  energy  used  for  the  two-ATF  model,  equation  (133),  is 


used  again: 

f  =  1/2E6?  +  l/2ai^i2  + 

~  “  67(^7  ~  Ct3^i^2  (166) 

Equation  (134)  is  the  equation  of  the  stress  field: 

a  =  df/d^=  E6  -  <5,^1  +  67^7  (167) 

The  affinities  are  given  by  equations  (135)  and  (136): 
A,  =  -^f/a^  =  <5i€  -  +  a, ^2  (168) 

Aj  =  ~  67^  ~  Ct2^2  0*3^1  (169) 

Balancing  the  linear  momentum  equation  gives 

dcr/ax  =  pu  =  Eu"  -  +  67(7'  (170) 

Recall  from  equation  (137)  that 

a^i/at  =  -Bi(^  -  }i)  i  =  1,2  (171) 

Allowing  fractional  derivatives  in  the  irreversible 
thermodynamic  equilibrium  equations  of  the  time  rate  of 


change  of  ^1  and  equations  (141)  and  (142),  results  in 
the  following: 

a'^^i/a't'^  =  ~Bi^i  +  Bi5i€/ai  Bio<3^z/ai  (172) 

a'^^?/a‘t'‘  ~  ~B?^2  +  B?(52^/a?  +  B2Ct3^i/a2  (173) 

Equations  (172)  and  (173)  together  with  equation  (170)  are 

the  constitutive  equations  for  a  uniaxial  rod  modeled  with 
two  augmenting  thermodynamic  fields  with  general  order 
derivatives.  These  equations  are  three  partial  differential 
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equations,  coupled  in  u,  and  ^2,  and  can  be  solved  more 

easily  in  the  Fourier  domain. 


Damping  Ratio.  Taking  the  Fourier  transform  of 
equations  (170),  (172),  and  (173),  with  iu)  introduced 
through  the  time  derivatives  and  ik  introduced  through 
derivatives  with  respect  to  x,  gives  three  equations  which 
may  be  written  in  matrix  form  as 


-pu  ^  +  Ek  2 

ik(5i 

ik(52 

u 

iB 1 k  d 1 

-Oi  tBi  +  (i(x>)q] 

BjOtj 

< 

Zi 

iB2k(52 

B20t3 

-02  (B2  +  (ica)>'] 

Z2 

The  determinant  of  the  above  matrix  must  equal  zero  in  order 
to  solve  the  linear  system  of  equations  which  it  represents. 
The  determinant  of  equation  (174)  is  then 


(-pa>2+Ek2 )  {  (aia  2 )  (Bi +  ( icj)  q  ]  [B2  +  (  iw) '■  ]  “  816203^) 

-  (ik(5i  ){(-a2)  [iBiW  1]  [B2  +  (ia;)^)  -  iBiB2ka3(j2> 

+  (ik(52)<BiB2ka3(5i  -  (-oti )  [iB2k(52]  [Bi  +  (iw)<j]  >=0  (175) 


After  multiplication  of  terms  inside  the  {  >  brackets, 
equation  (175)  can  be  written  as 

(-P(j^+Ek2 )  { (a  10  2 )  [Bi B7+B1  ( iu>)  'J+B7  ( ia^)  r  +  ( ia^)  q  ♦  r  ] 

~  818203^)  —  (k^  <5 1 )  [BiB20*2(J  1 +Bi  ( iw) '"O  2  (5 1 

+BiB203(J2]  ”  (k^i 2 )  IB1B2O3  (Ji +81820*1  (5  2 

+B2  (icj)qoi  (52]  =  0  (176) 

Before  solving  this  set  of  equations,  the  terms 
involving  (iuj)q  and  (iu>)r  must  be  written  in  a  form  that 
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allows  for  easier  algebraic  manipulation.  Recall  from  the 
appendix  that  i*?  =  cos(q7r/2)  i  sin(q7r/2)  ,  and  allow 
q  +  r  =  a 

q7r/2  =  X 

ttt/I  -  y 

S7r/2  =  z  (177) 

Equation  (176)  may  now  be  written  in  the  following  form: 

(-P(j2+Ek2 )  {aia2  [BiB2+Bi<jrcos(y)+B2(^<icos(x)+tL»»cos(z) 

+  i  (Bitj'^sin(y)+B2tJ'Jsin(x)+£j8sin(  z ) )  ]  -BiB20<3  3} 

~  k 2  (B 1 B 2O 3  b\b  1 +B 1 B 20^2  5 1  ^ ■*‘B i(j'‘‘0<25i  ^cos(y) 

+  iBitjrot^a  1 2sin(y )  ]  -  k2  [BiB20»3  (5i  5  2+BiB2ai  (5  2  2 

+  B26j'50<i  (J2  ?cos  (x)-iB2tj9at  d 2  ^sin(x)  ]  =  0  (178) 

Dividing  this  equation  by  the  term  which  multiplies 
(-pu)^  +  Ek2)  results  in  an  equation  of  the  form 

-pu;2  +  Ek2  -  k2(Ni/Di)  =  0  (179) 

where 

Ni  =  2BiB20(3()  1()  2  +  BiB20*23i2  +  BlB2aid2^ 

+  Bi<*2ra2  2co3(y)  +  B2a»‘2ai d  2  ^ cos (x) 

+  i  [Bia;ra2  <5i  ^sin(y )  +  B2tJ*’ai5  2  ^sin(x)  ]  (180) 

and 

Di  *  BiB2(aia2  -032)  +  aia2  (Biw^ cos (y )  +  B2a>‘Jcos(x) 

+  tJ*cos(z)]  +  iaia2  [BiCi;>^sin(y)  +  B2Ci;'isln(x) 

+  Ci2»sin(x)]  (101) 

Factoring  EB1B2OI1CH2  from  the  numerator,  Ni,  and 
recalling  the  definitions  of  ^1,  A?,  and  ds  from  equation 
(150)  results  in 
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Ni  »  EBiB20fia2{43  4i[l  +  (tJ>"/B2  )co5(y)] 

+  A7[1  +  (o^-J/Bi  )cos(x)l  +  i  [(4itJ''/B2  )sin(y) 

+  (426l»^/Bi  )sin(x) ]  >  (182) 

Factoring  BiB20t:a2  from  the  denominator,  Di ,  and 
recalling  that  X  =  1  -  a3  2/(oiiO»2)  gives 

Di  =  BiB2aiOt2{y+  (tj<J/Bi  )cos(x)  +  (tj^ /B2  )cos(y ) 

+  (cj*/BiB2  )cos(z)  +  i  ((o^-J/Bi  )ain(x)  +  (tjr /B2  )8in(y) 

+  (o^«/BiB2  )ain(z)]  >  (183) 

Dividing  equation  (179)  by  p  and  rearranging,  with  the 
EBiB2C(iO(2  canceling,  leaves  an  equation  from  which  the 
complex  modulus  can  be  extracted: 

-u>2  +  E(k2/P)[l  -  N2/D2]  =  0  (184) 

where  N2  and  D2  are  defined  by  equations  (182)  and  (183) 
divided  by  BiB20fioe2,  repectively.  The  equation  inside  the 
(  ]  brackets  is  the  complex  modulus  divided  by  E,  and  the 
loss  factor  can  be  derived  from  equation  (4).  To  do  this, 
however,  both  N2  and  D2  must  be  multiplied  by  the  complex 
conjugate  of  D2 : 

D2  =  y  +  (w<i/Bi  )cos(x)  +  (oi"^ /B2  )cos (y ) 

+  (w»/BiB2  )cos(z)  -  i  ( (w'J/Bi  )sin(x) 

+  (wr/B2)8in(y)  +  (a;«/BiB2  )sin(z)  ]  (185) 

In  order  to  simplify  the  resulting  equations  throe 
trigonometric  relations  were  required: 

sin(x  +  y)  =  sin(x)cos(y)  +  cos(x)oin(y) 
cos(x  +  y)  =  cos(x)cos(y)  +  sin(x)sin(y) 

A  =  A  sin2(x)  +  A  cos2(x)  (186) 
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Following  extensive  algebraic  manipulation  and 
employing  the  above  trigonometric  relations  results  in  the 
following  expressions  for  the  numerator,  N2,  and  the 
denominator,  D2,  of  equation  (184): 

Ni  (1/E)  -  Ax  [y+  (wr/B2)2  +  (a><J/Bi)cos(x) 

♦  (  y  l)(02''/B2)C0S(y)  +  (<J*/BiB2  )cos(z) 

+  (<j*/BiB2  )coa(x-y)  +  (a>(  r  > /B1B2  ^  )cos{z-y)] 

+  42  [y+  (cj<J/Bi)2  +  (y+  1)  (cjq/Bi  )cos(x) 

+  ((J^'/B?  )cos(y )  +  (tj*/BiB2  )cos(z) 

+  (cl>»/BiB2  )cos(x-y)  +  (a>' <!♦»  > /Bi  ^B2  )cos(z-x)] 

+  43  tv  +  (cj<J/Bi  )cos(x)  +  ((j'‘/B2  )cos(y) 

+  (tj»/BiB2  )cos(z)] 

+  i  {  4i  [-(u>VBi)sin(x)  +  (y  -  l)(wr/B2)sin(y) 

-  ((j»/BiB2  )sin(z)  +  (a;* /BiB:  )sin(y-x) 

+  (cj<  > /B1B2  nsin(y-z)] 

+  42  [O'  -  1)  (tj<i/Bi  )sin(x)  -  (02*- /B2  )sin(y) 

-  (cl»*/BiB2  )sin(z)  +  (a»«/BiB2  )sin(x-y ) 

+  ((j(  ) /Bi  2B2  )sin(x-z)] 

-43  ((w<J/Bi  )sin(x)  +  (t.;r/B2  )sin(y) 

+  (a2«/BiB2)sin(z)]}  (187) 

D2  «  yz  +  2y  ((w<J/Bi  )cos(x)  +  (tJ*'/B2  )cos(y) 

+  (a2*/BiB2  )cos(z)]  +  (cj<J/Bi)2  +  +  (tJ»/BiB2)2 

+  (2a2*/BiB2)  [cos(x-y)  +  (a;9/Bi  )cos(z-x) 

+  (02^/32 ) cos (y-z)]  (188) 
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The  loss  factor  becomes 


rj  »  Im[E*(w)]/J?otE*(w)] 

=  (Im[D2-N2]/D2)/(J?e[D2-N2]/D2)  (109) 

Again,  the  assumptions  of  small  damping  and  Ai  <<  1 
A2  <<  1  .  and  As  <<  1  .  loading  to  equation  (156)  are 

made  hero.  Two  other  relationships  help  further  reduce  the 
resulting  loss  factor.  Recalling  the  definitions  of  x,  y, 
and  z  given  in  equation  (177),  gives  the  following 
relationships : 

z  -  y  =  (q+r)(7r/2)  -  rn/l  =  q7r/2  =  x 

z  -  X  =  (q+r)(7r/2)  -  q7r/2  =  ttt/I  =  y  (190) 

Using  trigonometric  relationships  for  the  sine  and  cosine 
functions  results  in 

sin(x-z)  *  -  sin(z-x)  =  -  sin(y) 
sin(y-z)  =  -  sin(z-y)  =  -  sin(x) 
cos(x-z)  =  cos(z-x)  =  cos(y) 

cos(y-z)  =  cos(z-y)  =  cos(x)  (191) 

The  loss  factor  then  becomes 

V  =  Nlosa/Dloss  (192) 

where 


Nioi.  {((w’/Bi)  +  (a;(  *■♦•  ) /B1B2  n  ]  9in(x) 

+  (cj«/BiB2  )  [sin(z)  +  sin(x-y)] 

+  (1  -  y)(a;r/B2)sin(y)>  +  ^2  {[(1  -  y  )  (a><i /Bi  )sin(x) 
+  [((jr/B2)  +  (<^<  <I  ) /Bi  2B2  )]  sin(y)  , 

+  ((j« /B1B2 )  [sin(z)  +  sin(y-x)]> 

+  [  (tjq/Bi  )sin(x)  +  (a;r/B2  )oin(y ) 

(a;*/BiB2  )sin(z) ) 
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(193) 


and 


•  Dioe.  =  y  2  +  2v[(a>i/Bi  )cos(x)  +  (tjr/Bz  )cos(y) 

+  (tj»/BiB2)c08(z)]  +  (a;<i/Bi)2  +  (uVB2)2 
+  (a;»/BiB2)2  +  (2t*>«/BiB2 )  [cos(x-y)  +  (tJ'J/Bi  )cos(y) 

•  +  (a2r/B2)coa(x)]  (194) 

With  the  assumption  of  small  damping,  the  damping  ratio  is 
equal  to  half  the  loss  factor. 

•  As  an  added  check  on  the  algebra,  the  loss  factor  of 
the  two-ATF  model  with  fractional  derivatives,  equation 
(192),  should  reduce  to  the  loss  factor  of  the  single-ATF 

•  model  with  fractional  derivatives,  twice  equation  (116),  in 
the  same  manner  equation  (156)  reduced  to  equation  (41)  when 
the  two-ATF  model  without  fractional  derivatives  was 

•  compared  to  the  single-ATF  model.  Setting  Qts  =  (52  =  0  , 
which  in  turn  drives  A2  =  da  =  0  and  y  =  1  ,  in  equation 
(192),  then  multiplying  the  numerator  and  denominator  by  B2^ 

•  results  in  the  following: 

^  =  Nck/Dck  (195) 

where 

•  Nek  =  di  {Ba  z  (<J<*/Bi  )sin(x) 

+  B2  (cj*/Bi )  (8in(z)  +  sin(x-y)  -f  u>rsin(x)  1 }  (196) 

and 

•  Dck  «  B22  (1  +  2((j<i/Bi  )co8(x)  +  ((^i/B:)2] 

B2  [2tj*‘C08(y)  +  2(cj»/Bi  )co8(z)  2(a>»/Bi  )co8(x-y) 

2(tj(  q**  ) /Bi  2  )co8(y)l  +  2(02*  > /Bi  )cos(x) 

•  (o^O^  +  (w»/Bi  )^  (197) 
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Taking  the  limit  of  equation  (195)  as  B?  goes  to  infinity 
results  in  an  indeterminant  form.  L'Hospital's  rule  is 
invoked  twice,  and  the  resulting  limit  gives 

)sin(x) 

7]  s  -  (198) 

1  +  (2tjq/Bi  )cos(x)  +  (tj‘i/Bi)2 
which  exactly  equals  twice  equation  (114),  and  the  loss 
factor,  therefore,  reduces  to  that  of  the  single-ATF 
fractional  derivative  model. 

Characteristics  of  the  Two-ATF  Model  With  Fractional 

Derivatives .  Figure  20  is  a  plot  of  the  loss  factor  versus 

frequency  for  the  two-ATF  model  with  the  values  of  the 

parameters  as  follows: 

q  =  0.8  Bi  =  1  Ai  =  0.1 

r  *  0.4  B2  =50  A?  =  0.05 

A3  “  5  X  10-2 
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Figure  20.  Loss  Factor  vs.  Frequency  for  the  Two-ATF  Model 

with  Fractional  Derivatives 


Figure  20  shows  a  similar  "double  hump"  characteristic 
to  that  of  the  two-ATF  model  without  fractional  derivatives. 
The  first  maximum  falls  under  the  value  of  frequency 
corresponding  to  the  value  of  Bi .  However,  unlike  the 
two-ATF  model  without  fractional  derivatives,  the  second 
maximum  would  appear  not  to  correspond  to  B2.  Comparing 
equation  (192)  to  equation  (156)  shows  that  the  Bi  of 
equation  (156)  are  raised  to  the  same  power  as  cj.  In 
equation  (192),  however,  this  is  not  the  case.  To  find  the 
value  of  frequency  at  which  the  maxima  for  the  loss  factor 
corresponding  the  two-ATF  model  with  fractional  derivatives 
will  occur,  the  values  of  the  Bi  must  be  raised  to  the  power 
of  inverse  of  the  respective  fraction  to  which  the 
corresponding  d  is  raised.  For  example, 

Bi  =  1 ,  q  =  0 . 8 

==>(j=(l)<i/o.e)  =  (i)i.25  =  i  (199) 

which  is  the  value  of  frequency  corresponding  to  the  first 
maximum.  Likewise,  for 

B?  =  50,  r  =  0.4 

==>  (jJ  =  (50)(i/o.4)  =  (50)2-5  =  17670  (200' 

Again,  the  value  corresponds  to  the  frequency  of  the  second 
maximum  in  figure  20. 

Notice  also  that  the  slopes  of  the  two  portions  of  the 
curve  are  different,  with  the  slopes  of  the  first  maximum 
being  greater  than  those  of  the  second  maximum. 


Equations  (166),  (172),  and  (173)  are  symmetric  with 
respect  to  the  two  ATFs,  as  well  as  the  other  material 
constants.  This  fact  leads  to  the  conclusion  that  the  loss 
factor  must  also  be  symmetric  with  respect  to  the  parameters 
of  the  model.  Exchanging  the  values  of  all  parameters  with 
subscript  1  for  those  with  subscript  2,  e.g.,  Bi  =  B?  and 
B2  =  Bi  ,  as  well  as  exchanging  the  values  of  the 
fractional  powers,  q  and  r,  results  in  a  plot  of  the  loss 
factor  that  is  exactly  the  same  as  figure  20.  This,  then, 
is  an  added  check  on  the  algebra  performed  to  arrive  at 
equation  ( 192 ) . 

Figures  21  through  34  show  the  loss  factor  versus 
frequency  as  various  parameters  of  the  two-ATF  model  with 
fractional  derivatives  are  varied. 

The  value  of  q  is  reduced  to  0.1  in  figure  21.  The 
first  maximum  is  no  longer  visible,  and  the  value  of  the 
second  maximum  has  increased.  When  q  is  increased  to  0.95 
in  figure  22,  the  value  of  the  peak  is  increased,  while  the 
maximum  of  the  second  peak  remains  the  same.  The  slope  of 
the  first  "hump"  or  curve  is  clearly  greater  as  well,  as  is 
evidenced  by  the  narrowing  of  the  peak,  and  the  deepening  of 
the  minimum  as  the  first  curve  leads  into  the  second  curve. 

Values  of  r  are  varied  in  figures  23  and  24.  First,  r 
is  reduced  to  0.1,  and  the  second  peak  shifts  to  the  right. 
The  maximum  value  of  the  first  peak  remains  unchanged.  As  r 
is  increased  to  0.65,  the  second  peak  shifts  towards  the 


Figure  22.  Fractional  Derivative  2-ATF  Model  with  q=0.95 
{r=0.4,  B]=l,  B?*50,  Ai=0.1,  d?=0,05,  A3*'/5  X  IQ-?) 


Loss  Factor 
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first  peak  on  the  left,  and  the  value  of  the  second  maximum 
increases.  Once  again,  the  value  of  the  first  maximum 
remains  constant.  It  will  also  bo  noted  that  the  value  of 
the  slope  of  the  second  curve  has  increased  as  well. 

Figure  25  shows  the  plot  of  the  loss  factor  versus 
frequency  when  Bi  is  decreased  to  0.3.  The  peak  now 
corresponds  to  a  frequency  of  0.22  (u>  =  (0.3)i-25),  with  the 
slopes  and  maximum  the  same  as  originally  plotted  —  the 
first  maximum  has  just  shifted  to  the  left.  The  second  peak 


Figure  25.  Fractional  Derivative  2-ATF  Model  with  Bi=0.3 
(q=0.8,  r=0.4,  B2=50.  Ai=0.1,  A7*0.05,  X  lO'Z) 

is  now  also  more  clearly  discernable.  When  Bi  is  Increased 
to  30  in  figure  26  the  curve  shifts  to  the  right  with  the 
maximum  corresponding  to  a  frequency  of  70.  The  second 
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Figure  26.  Fractional  Derivative  2-ATF  Model  with  Bi=30 
(q=0.8,  r=0.4,  B2  =  50,  A?*0.05,  A3  =  v^  X  lO'Z) 


maximum  is  now  almost  a  part  of  the  first  curve,  but  can 
still  be  seen  as  an  inflection  point. 

With  B2  reduced  to  5  in  figure  27  the  plot  of  loss 
factor  versus  frequency  shows  the  shift  of  the  second 
maximum  to  the  left,  with  the  maximum  now  an  inflection 
point  on  the  curve  at  a  frequency  of  56.  The  slope 
following  the  inflection  point  clearly  differs  from  the 
elope  preceeding  the  point,  showing  the  elope  dependence  or. 
the  fractional  powers  of  the  model.  Likewise,  as  the  value 
of  B2  is  increased  to  1000  in  figure  26,  the  maximum  of  the 
second  curve  shifts  to  the  right.  The  maximum  should 
correspond  to  a  frequency  of  3  x  10^ . 
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Figure  27.  Fractional  Derivative  2-ATF  Model  with  B?=5 
(g=0.8,  r=0.4,  Bi=l,  Ai=0.1,  A?*0.05,  ^3-\f5  X  10-z) 


^  Figure  28.  Fractional  Derivative  2-ATF  Model  with  B2“1000 

(q*0.8,  r=0.4,  Bi=l,  di«0.1,  A2«0.05,  X  10-2) 


'3 


The  value  of  affects  the  value  of  the  maximum  of  the 
first  peak,  as  can  be  seen  in  figures  29  and  30.  Reducing 
to  0.01,  and  order  of  magnitude,  reduces  the  maximum  by 
almost  an  order  of  magnitude  in  figure  29.  Increasing  the 
value  of  has  the  opposite  effect  of  increasing  the 
magnitude  of  the  peak  in  figure  30.  In  both  cases,  however, 
the  value  of  the  second  maxim\jun  remains  unaffected. 


Figure  29.  Fractional  Derivative  2-ATF  Model  with  A] =0.01 
(q=0.8,  r=0.4,  Bi  =  l,  B2=50,  A2=0.05,  A3  =  n/5  X  lO"?) 

Figures  31  and  32  show  the  effects  of  varying  a?.  If 
A;  is  decreased  as  in  figure  31,  the  value  of  the  second 
peak's  magnitude  decreases  accordingly.  Increasing  A?  to 
0.5  in  figure  32  increases  not  only  the  magnitude  of  the 
second  peak  quite  drastically,  but  also  that  of  the  first 
peak,  although  to  a  much  weaker  extent. 
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Figure  30.  Fractional  Derivative  2-ATF  Model  with  Ai=0.9 
(qsO.8,  r=0.4,  Bi=l,  B?=50,  A7=0.05,  X  10-2) 
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Figure  32.  Fractional  Derivative  2-ATF  Model  with  Aj=0.5 
(q=0.8,  r=0.4,  Bi=l,  B7=50,  Ai=0.1,  A3  =  n/5  X  10-2) 


Finally,  the  value  of  A3  is  varied  in  figures  33  and 
34.  Reducing  A3  toVoTs  X  10-3  by  reducing  C’  an  order  of 
magnitude  does  little  to  the  overall  curve  except  for 
reducing  the  value  of  the  magnitude  of  the  first  peak  only 
slightly.  Increasing  Ai  to  \/40 . 5  X  10"  2  by  increasing  C*  to 
0.9,  however,  increases  the  value  of  the  first  peak  by  an 
order  of  magnitude,  and  shifts  it  to  the  left.  The  second 
peak,  in  turn,  is  only  affected  slightly,  with  a  small 
increase  in  the  magnitude  visible. 

The  effects  of  variations  in  the  model  parameters  are 
summarized  in  figure  35.  Increasing  the  parameter  will 
shift  the  curve  in  the  direction  shown. 
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Figure  33.  Two-ATF  Fractional  Derivative  Model 
with  d3=v^^  X  10-3 

(a=0.8,  r=0.4,  Bi=l,  B?*50.  di=0.1,  ^7=0. 05) 
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Figure  34.  Two-ATF  Fractional  Derivative  Model 
with  43  =  v/40.5  X  10-2 

(q=0.8,  r=0.4,  Bi*l,  B7=50,  d]=0.1,  47*0.05) 
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Figure  35.  Effect  of  Parameter  Variation  on  Two-ATF 
Model  with  Fractional  Derivatives 
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In  this  study,  the  development  of  a  model  of  material 
damping  using  augmenting  thermodynamic  fields  (ATF) ,  wherein 
the  equations  of  thermodynamic  equilibrium  for  these  ATFs 
had  derivatives  of  integer  order,  was  reviewed.  This  model 
was  originally  introduced  by  Lesieutre,  formerly  of  UCLA. 

The  model  was  shown  to  be  able  to  reproduce  the  material 
damping  of  classical  thermoelastic  and  viscoelastic  theory, 
with  no  assumptions  made  as  to  the  magnitude  of  the  the 
parameter  A,  given  that  E  is  taken  to  be  the  adiabatic 
modulus  of  elasticity  for  the  thermoelastic  case,  or  the 
glassy  modulus  for  the  viscoelastic  case. 

As  an  aside,  if  A  is  small  as  Lesieutre  assumed,  then 
the  modulus  of  elasticity  in  the  ATF  model  formulation,  E, 
is  approximately  equal  to  the  isothermal  modulus  of 
elasticity,  Et,  of  the  thermoelastic  model,  or  equal  to  Eo, 
the  Young's  modulus  of  the  viscoelastic  model. 

The  model  of  material  damping  was  also  developed 
wherein  the  thermodynamic  equilibrium  equations  of  the 
augmenting  thermodynamic  fields  were  allowed  to  have 
derivatives  of  fractional  order.  The  ATF  model  with 
fractional  derivatives  was  then  compared  to  Bagley  and 
Torvik's  4-parameter  viscoelastic  model,  in  which  fractional 
derivatives  are  allowed  as  well.  Comparison  was  a  one  to 
one  match,  given  the  fractions  of  each  model's  derivatives 
were  equal.  This  comparison  also  showed  that  the  A 
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parameter  of  the  fractional  order  ATF  model  was  less  than 
unity,  Moreover,  in  many  materials  A  approached  unity, 
i.e.,  A  >  0.999  ,  since  the  glassy  modulus  for  many 

materials  is  three  orders  of  magnitude  greater  than  the 
rubbery  modulus.  Owing  to  the  fact  that  Bagley  and  Torvik 
have  used  the  4-parameter  model  to  accurately  model  many 
materials  over  broad  frequency  ranges,  assumptions  that 
d  <<  1  ,  as  Lesieutre  made  in  his  original  model,  cannot  be 

made  when  using  the  ATF  model  to  describe  viscoelastic 
behavior  in  material  damping. 

The  ATF  model  and  the  ATF  model  with  fractional  order 
derivatives,  a  more  general  case  of  the  ATF  model,  was  shown 
to  lead  to  expressions  for  the  loss  factor  and  damping  ratio 
that  are  strong  functions  of  frequency,  for  the  case  of  the 
uniaxial  rod.  This  was  the  desired  result.  More 
importantly,  the  ATF  model  accurately  describes  actual 
behavior  of  many  materials,  because  of  the  correspondence 
between  it  and  the  models  of  thermoelasticity  and 
viscoelasticity. 

In  all  comparisons,  the  modulus  of  elasticity,  E,  of 
the  ATF  model  differed  from  the  moduli  of  the  respective 
models,  Et  and  Eo ,  by  the  factor  (1-4).  In  the  case  of 
viscoelasticity,  the  modulus  E  is  found  to  be  the  glassy 
modulus,  Ea.  with  the  relationship  between  the  rubbery 
modulus,  Eo,  and  the  glassy  modulus  being  that  of 
Eo  =  (1  -  4)Ea  A  similar  relation  was  shown  to  exist 


between  the  isothermal  and  adiabatic  moduli  of  the 


thermoelastic  model,  i.e.,  Ea .  the  adiabatic  modulus  of 
elasticity  is  related  to  Et .  the  isothermal  modulus  of 
elasticity  by  the  equation  Ea  =  (1  -4)Et. 

Recommendations  for  further  study  include  examination 
of  the  linearity  effects  of  this  model.  Determination  of 
whether  the  sum  of  results  obtained  from  a  fractional  order 
ATF  modeling  of  viscoelasticity  will  describe  material 
damping  of  materials  in  which  both  thermoelasticity  and 
viscoelasticity  are  present  need  be  made.  It  would  appear 
to  be  so,  since  the  fractional  order  ATF  model  is  linear. 
This  would  also  include  experimental  determination  of  the  a  3 
coupling  parameter  between  strain  and  temperature. 

It  is  recommended  that  the  effects  of  higher  order 
derivatives  in  the  diffusion  laws  on  the  constitutive  laws 
be  determined.  This  will  lead  to  an  understanding  of  the 
exact  relationship  of  the  fractional  order  ATF  model  to 
thermoelasticity  and  viscoelasticity.  The  expansion  of  the 
model  to  three  dimensions  should  present  no  particular 
difficulty. 

Possible  future  research  could  also  include  an 
examination  of  energy  methods  as  a  means  of  deriving  the  ATF 
model . 

Finally,  it  is  recommended  that  numerical  methods  for 
the  solution  of  the  fractional  order  ATF  model  be  ananlyzed. 
Application  of  numerical  analysis  to  the  fractional  order 
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ATF  model  will  allow  for 
model  in  the  description 


further  ease  in  the  use  of  the 
of  material  damping. 
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Appendix  A: 

Complex  Numbers  (Kreyszig,  1988:721-729) 


Define  a  complex  number,  z,  as 
z  =  X  +  iy 

where  both  x  and  y  are  real  numbers,  and 

i2  =  -1 

is  the  definition  of  an  "imaginary  number." 

The  complex  conjugate,  z,  of  z  is  defined  as 
z  =  X  -  iy 

From  these  definitions,  several  relationships  are 
developed: 

Re  z  -  -x.  -  l/2(z  +  z) 

Im  z  =  y  =  [l/(2i)](z  -  z) 


(A.l) 


(A. 2) 


(A. 3) 


(A. 4^ 

(A. 5) 


Zl\  Zi 


(A. 6) 


22/  22 


The  polar  representation  of  the  complex  number  z  is 
often  practical  to  use,  and  easier  to  manipulate 
mathematically.  Define 
X  =  r  cos(^) 


y  =  r  cos(^) 


(A. 7) 


By  solving  these  equations  simultaneously,  it  can  be  shown 


that 


r  =  (x2  +  y2)i^‘ 
0  =  arctan(y/x) 


(A. 8) 


Substituting  equations  (A. 7)  into  equation  (A.l)  gives 


2  =  r(cos(0)  +  i  sin(^)] 


(A. 9) 


83 


It  can  also  be  shown  mathematically  that 
zn  =  rr'[cos(n(y)  +  i  sin(n^)] 


(A. 10) 


One  of  the  most  important  complex  analytical  functions 
is  the  exponential  function,  defined  as 

ez  =  e5<[cos(y)  +  i  sin{y)]  (A. 11) 

Substituting  z  =  iy  into  equation  (A. 11)  gives  the 
so-called  Euler  Formula: 

eiy  =  cos(y)  +  i  sin(y)  (A. 12) 

Comparing  the  right  hand  side  of  equation  (A. 9)  to  that 
of  equation  (A. 12)  gives 

eiO  =  cQsiO)  +  i  sin(^)  (A. 13) 

and  together  with  equation  (A. 9)  gives  the  polar  form  of  the 
exponential  function: 

2  =  r  ei«  (A.  14) 

Solving  equation  (A. 14)  when  z  =  i  gives 
i  =  r[cos(^)  +  i  sin(^)]  (A. 15) 

For  the  equality  of  equation  (A. 15)  to  hold  the  following 
must  be  true: 

r  =  (X2  +  y2)l/^  =  (02  +  12)1/2  =  1 

sin(^)  =  1  (A. 16) 

Therefore , 

6  =  (n7r)/2,  n=l ,  2,3.  .  .  (A.  17) 

So,  recalling  equation  (A. 13), 

i  =  cos(7r/2)  +  i  sin(7r/2)  =  ei('^/2)  (A. 18) 

for  the  principal  value  (- tt  <  $  <7T). 
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For 
z  = 


z  =  ( i ) f  where  0  <  f  <  1 
if  =  (ei(rr/2))f  =  (A. 19) 


From  equations  (A. 10)  and  (A. 18) 

if  =  cos(^)  +  i  sin(0)  =  cos(f7r/2)  +  i  sin(f7r/2)  (A. 20) 


Appendix  B: 

Specific  Heat  Comparisons  (Torvik,  1990) 


Let  the  free  energy  of  a  solid  undergoing  a 
one-dimensional  thermoelastic  deformation  be 

F  =  E  -  TS  (B.l) 

where 

E  =  stored  elastic  energy 
T  =  absolute  temperature 
S  =  entropy. 

By  the  combined  first  and  second  laws  of  thermodynamics  for 
this  deformation 

dE  =  TdS  +  a (B.  2) 

But 

dF  =  dE  -  TdS  -  SdT  (B.3) 

So 

dF  =  Ode  -  SdT  (B.4) 

Thus 

aF/aT)^  =  -S  (B.5) 

and 

dF/de)T  =  o  (B.6) 

Hence  it  follows  that 

as/aeji  =  -a2F/(a€aT)  =  -a<T/aT),  (b.?) 

Recall  that,  for  constant  pressure, 

CpdT  =  dQ  =  TdS  (B.8) 

or 

Cp  =  T  ds/dr)^  (B.9) 

and,  for  constant  volume  (normal  strain). 


CvAT  =  AQ  =  TAS  (B.  10) 

or 

Cv  =  T  ds/dj),  (B.ll) 

Now  consider  an  incremental  change  in  strain,  at 
constant  entropy: 

do  =  ^o/^€)^  d€  +  do/dT)t  dT  (B.12) 

So 

do/de  )  s  =  do/d€)T  +  dT/de)s  do/dT^^  (B.13) 

Recall  also  that 

ax/0y)j  ay/0z)>  az/ax)y  =  -i  (b.i4) 

Hence 

dT/de)s  di/dS  ]  T  as/ar),  =  -1  (B.15) 

So 

dT/dc)%  =  -  ds/di)i  ai/as),  (B.16) 

Also, 

^o/^€)^  de/d’l)„  d1/do)^  =  -l  (B.17) 

So 

acr/ai)*  =  -  do/d€)T  d€/dT:)o  (B.18) 

Substituting  equation  (B.16)  into  equation  (B.13),  then 
using  equations  (B.7),  (B.8)  and  (B.ll)  we  find 

a<7/ae)s  =  do/df)?  +  td€/do)T]HdT/d€)}^  T/Cv  (B.19) 

or 

do/de)^  =^o/^€)^  +  [aa/ae)!]  ^tae/ax)^]  2  T/Cv  (B.20) 
or 

E»  =  Et  +  Et^O^T/Cv  (B.21) 

where 
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Ea  is  the  adiabatic  Young's  modulus, 

Et  is  the  isothermal  modulus,  and 

ot  is  the  linear  coefficient  of  thermal  expansion. 


Thus , 


Ea  =  Et[1  +  (EtO^TI/CvI 

Note  that  Cv  here  is  the  energy  change,  per  unit 


volume,  for  a  unit  temperature  change. 
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ABSTRACT 


This  study  developed  a  material  damping  model  using  the 
concept  of  augmenting  thermodynamic  fields  |aTFJ-t  wherein  the 
equations  of  thermodynamic  equilibrium  are  allowed  to  have 
derivatives  of  fractional  order. j  This  effort  seeks  to  expand 
the  applicability  of  the  ATF  model  of  Lesieutrei  formerly  of 
UCLAt  in  which  ATFs  interact  with  the  mechanical  displacement 
fields • 

Current  models  of  material  damping  cannot  predict  well 
the  dependency  of  damping  on  frequency-  Two  newer  models-^ 
capable  of  predicting  this  frequency  dependence -I'^i-are  discussed- 
They  include  the  ATF  model-i  and  Bagley  and  Torvik's-i  AFIT-. 
4-parameter  modeli  which  allows  fractional  derivatives  in  the 
description  of  viscoelastic  materials- 

This  research  effort  applies  fractional-order  derivatives 
to  the  ATF  model-  Coupled  material  constitutive  relations 
are  developed  using  the  concept  of  augmenting  thermodynamic 
fieldST  with  non^integer  differentials  allowed  in  the  resulting 
partial  differential  equations. 

The  complex  modulus  that  results  from  solution  of  these 
partial  differential  equations  is  compared  to  the  complex 
moduli  of  thermoelasticity -I  integerf  order  viscoelasticity-) 
and  viscoelasticity  with  fractional  derivatives  -Cthe 
4-parameter  model!  for  the  case  of  a  uniaxial  rod-  In  each 
casei  the  fractional -order  ATF  model  reduced  to  the  respective 
modeli  and  therefore-i^accurately  describes  the  damping 
mechanism  resulting  from  each  of  these  models-  > 


